{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PHydro-cover-small.png\">\n",
    "*This is the Jupyter notebook version of the [Python in Hydrology](http://www.greenteapress.com/pythonhydro/pythonhydro.html) by Sat Kumar Tomer.*\n",
    "*Source code is available at [code.google.com](https://code.google.com/archive/p/python-in-hydrology/source).*\n",
    "\n",
    "*The book is available under the [GNU Free Documentation License](http://www.gnu.org/copyleft/fdl.html). If you have comments, corrections or suggestions, please send email to satkumartomer@gmail.com.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Multivariate Distribution](10.02-Multivariate-Distribution.ipynb) | [Contents](Index.ipynb) | [Bias Correction](10.04-Bias-Correction.ipynb)>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10.3 Kriging\n",
    "\n",
    "克里金(Krigin)是水文学领域中广泛应用的一种插值方法。克里金优于其他插值方法，例如线性插值、径向函数法和样条插值法等，它提供了一种方法来估计插值参数，同时也提供了插值的不确定性。克里金方法包括两个步骤：(i)对变差图进行估计并对其进行理论变差图拟合，以及(ii)进行插值。经验变差图揭示了数据和采样间距的许多许多有趣的事情。研究这些方面是值得的。所以，首先我们开始通过不同类型的虚拟数据集来理解变差图。在本节中的例子来自P.K.Kitanidis的`Introduction to Geostatistics`一书。python库用于Kriginig是`ambhas.krige`。\n",
    "\n",
    "首先，我们考虑变异规模大于抽样规模的情况。考虑这个变量，\n",
    "\n",
    "<center>$z(x)=cos(2x/0.001)\\quad\\quad\\quad\\quad(10.1)$ </center>\n",
    "\n",
    "其中，$x$表示从范围在0到1之间的均分分布抽样100次。因为`ambhas.krige`库只适用于二维数据，我们将多生成一个除x之外的变量$y$。$y$值将是常量，只模拟一个维度的影响。生成虚拟数据的python代码是："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import required library\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.random.rand(100) # uniformly distributed x\n",
    "y = np.ones(x.shape)\n",
    "z = np.cos(2*np.pi*x/0.001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "让我们看看生成的虚拟数据。结果图如图10.3所示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X2wHXd93/H3R5JluIWAJV2IYvveaxNDMSQ16MaFMJMAxsU4M5YzJanJBeTUqQZct53QYbBHncI4aOIm03GemBoFHJtKAYPTFiWFOsbg0pYYfF2Mn1oj+UlW5eIbDDSMicHWt3/sXnR0dM65e87Z5/28ZnbOObt7zv52z+5+9/ewv1VEYGZmNq11VSfAzMzawQHFzMxy4YBiZma5cEAxM7NcOKCYmVkuHFDMzCwXDihmZpYLBxQzM8uFA4qZmeViQ9UJKNOWLVtiYWGh6mSYmTXKXXfd9dcRMbvWfJ0KKAsLCywvL1edDDOzRpH0WJb5XORlZma5cEAxM7NcOKCYmVkuHFDMzCwXlQYUSddLelLSfUOmS9IfSDoo6R5Jr+2ZtkPSgXTYUV6qzcxskKpzKDcAF4yY/jbgrHTYCfw7AEmbgA8Cfx84F/igpFMKTal1x759sLAA69Ylr/v2VZ0is0aoNKBExJeBp0bMsh34RCTuAF4saSvwVuDWiHgqIr4D3MrowJQfn2wm05Tttm8f7NwJjz0GEcnrzp31Ta9ZjVSdQ1nLqcDjPZ8Pp+OGjS9W3U42Pknnb9cuePrp48c9/XQy3sxGqntA0YBxMWL8iT8g7ZS0LGl5ZWVlutTU6WTjk3QxDh0ab3zRmnLRYEb9A8ph4PSez6cBR0aMP0FE7ImIxYhYnJ1ds+eA0ep0svFJuhhzc+ONL9Koi4ZxAo2DkpUlIiodgAXgviHTfgn4PEmO5HXA19Lxm4BHgFPS4RFg01rL2rZtW0xlfj4iObSPH+bnp/vdSUiD0yKVn5a11Gm7rWXv3oiZmePTOTOTjC/bsO22eXP2NNZpfayxgOXIcj7PMlNRA/BJ4AngRyS5jsuA9wDvSacL+AjwEHAvsNjz3X8MHEyHX8+yvKkDSp0OTp+ki7N3b7IdpeS1qnQOu2gYNgz676vcT+qyHW1qjQgoZQ9TB5SI+hwkPkm337BgMGwYlDutKifbtP3TRsoaUJTM2w2Li4vRqt6G9+1L6kwOHUrK+HfvhqWlqlNleVmtQ+mtK5uZgec/H7797RPnn5+HRx89ftzCQlL3Msj8fHH7zLDlDkqj1Z6kuyJica356l4pb6MsLSUH59GjyWtXg0lbK52XlmDPnuQkLCWve/bA7/9+Elh6zcwkwaHf7t0nzruqyJaBTWqIYfnJko1py5BLkZfVS1eLVsYpQlydd5y6l2mVXXeTZXu42HViuA7FAaUTmtQ4oWpl1qeUGeizLKurFx45yRpQXIdizbZuXXJ66CclRYF2TNn1GmXV8WVZL9fpTMV1KEVoa1l9k9XpRsS6G1SfMqzuJQ9l1fFlqa8pq06n4+cIB5SsmtTVSZeUfZKsqywnsmGV/E1vzJHloqKMCw+fI1yHkpnL6utrkgrqNlXMdr1+oC51KC0+R+BK+ZwDSpO6OrHB2nribfGJLLM6tPJq8Tkia0BxpXxWrtRrvrb+h26YUA9t3b9wpXz+XFbffG292c4NE+rB5wgHlMzaWqFZpbJbxLT1xOsTWT34HOE6FMso7/LnKuoz2lqHEtHOxgZWG7hS3gElN0WciKuqSPaJtzze1q2RNaC4Ut7WVkRloyuS221YT8ldKwJqCVfKW36KqMxua32GJZr0iGrLTaUBRdIFkh6UdFDSlQOmXyvp7nT4pqTv9kx7rmfa/nJT3jHDTvKbNk1eqe6K5HZra4s6G6mygCJpPcnjfd8GnA28Q9LZvfNExG9GxDkRcQ7wh8B/6Jn8g9VpEXFRaQnvokEn/5NOgr/5m8m7mXCLmOpN28pu1PedA62HsltSZqloKWIAXg/c0vP5KuCqEfN/BTi/5/P3x12mK+Wn0F/BunlzNZXqlo9pG1qs9f02t6gry7SNGnL8D6h7Ky/g7cDHej6/C/ijIfPOA08A63vGPQssA3cAF2dZ5kQBxS1VTrR37+Bg0pJuJjph2lZ2Wb7vY2dyeQSDHFtSNiGg/MqAgPKHQ+b9QP804KfS1zOBR4GXDfnuzjTwLM/NzY23FX2VdaJB26RJORSf5BLT9jvV4n6raiGPYJDjf5Q1oFRZKX8YOL3n82nAkSHzXgJ8sndERBxJXx8GbgdeM+iLEbEnIhYjYnF2dna8FLqlyokGbZNVda9Ud/fix0xbx9HGOpI6Pcskj0YNVfxHWaJOEQOwAXgYOAPYCHwDeNWA+V5BkgNRz7hTgJPT91uAA8DZay1z7CIvX4WdaNg2gfpf7btX3mOKrkNpmrqtTx77apfqUJI0ciHwTeAhYFc67mrgop55PgRc0/e9nwfuTYPQvcBlWZY3dkAZ50/tSlFKk0/KvkA4Xh6Vvm3Z5+u2X+cVDHL6jxoRUMoexg4oWf/Uul3dFKnJ61q3k4ZVp/9EOyzXXeXFRo0CtgNKHgElItuf2rUTVY129LE0ORhafgbtB8Nyr3kdw009ZlJZA4r78sqD+6Vqjn37koYFhw4llZO7d/tmyq4Z1jeddPxxnFffYy3o18x9eZWpjS1e2mppKenQ8ujR5LUhB7RllKWl1rCWUhHF9NzQodaiG6pOQCvs3j34CqTOTWjN2qY/J7DaLByODwxzc+U+qrdD/Zo5h5IH90tllqjyXo6sOYGyOybtUAmGA0peXJRiXVf1jaNZcwJlXwB2qGdtV8qbWT6KeBBbk5Y/SsMbg7hS3szKVXVdQZ1zAh0pwXBAMbN8VF1X4LrMyjmgmFk+6pBD6EhOoK4cUMwsH84hdJ7vQzGz/CwtOYB0mHMoZmaWCwcUMzPLhQOKmZnlwgHFzMxyUWlAkXSBpAclHZR05YDpl0pakXR3OvxGz7Qdkg6kw45yU25mZv0qa+UlaT3wEeB84DBwp6T9EfFA36w3RcQVfd/dBHwQWAQCuCv97ndKSLqZmQ1QZQ7lXOBgRDwcET8EPgVsz/jdtwK3RsRTaRC5FbigoHSamVkGVQaUU4HHez4fTsf1+4eS7pF0s6TTx/wuknZKWpa0vLKykke6rQpVdotuZplUGVA0YFx/18d/DixExM8CXwBuHOO7yciIPRGxGBGLs7OzEye2sco8ERe1rEHdor/rXXD55fn8ft05mFpDVBlQDgOn93w+DTjSO0NEfDsinkk//jGwLet3jXKfT1HksgY9OCkCrruu/SfXqp8xYjaGyp6HImkD8E3gPOD/AHcCvxYR9/fMszUinkjf/zLwgYh4XVopfxfw2nTW/wlsi4inRi2zc89DKfP5EEUua9265GQ6SB2edVGkOj/jwzoj6/NQKmvlFRHPSroCuAVYD1wfEfdLuhpYjoj9wD+XdBHwLPAUcGn63ack/RZJEAK4eq1g0kllPp+iyGUNewZ4Xr9fZ1U/Y8RsDJXehxIRn4uIl0fEyyJidzruX6fBhIi4KiJeFRF/LyLeFBH/u+e710fET6fDn1S1DpVaq2y9zOdTFLms3buT3muL+v06q/oZI2Zj8J3yTZWlbL3M51NMu6xRwXFpCd7znhODysaN8P3vt7uyuqj/0BX9VoSI6Mywbdu2aI35+YgklBw/zM8fP9/evck4KXndu7e4NE26rL17I2Zmjl+PmZkTv9/7+5s3R5x00trfaYO8/8Os29ssRVINseY5trJK+SrUplJ+376k5dKhQ0nRxe7d4z9DYlhFtZQ8ra5JJql4dmX15LztbExZK+Vd5FW2vJqBtqlsfZKKZ1dWT87bzgrigFK2QfdUPP10Mn4cdXh+d16GBcGI4eX7bQqoZfO2s4I4oJQtr6vDNj2/e1BwXDUsB9emgFo2bzsrSpaKlrYMtaiUz1qZ3jWrFc+Dts2w7VNmg4O28bZrt5z/X1wpf6JaVMqv1qH0FnvNzDQ3d5G3NjU2MKtCAecYV8rXVZuKqorg8n2z6eRVTzsBB5Ss8rwRbGkpaZ559Gjy6mByjMv3zaZTYSs+B5Qs3ONreZyDM5tOhbl816Fk4RvBzKwpXIdSc74RzMyaosJcfmXd1zfKsO7TXVFsZnW0tFRJMbFzKFm4otjMbE0OKFm4otjMbE2VBhRJF0h6UNJBSVcOmP4+SQ9IukfSbZLme6Y9J+nudNhfeGLd1NesGn52S2NUFlAkrQc+ArwNOBt4h6Sz+2b7OrAYET8L3Az8Ts+0H0TEOelwUSmJNrNyld1kf1DwckDLrLJmw5JeD3woIt6afr4KICJ+e8j8rwH+KCLekH7+fkS8YJxl1qLrFTPLrswm+4Oa227cmASyH/3o2LgOdpXUhGbDpwKP93w+nI4b5jLg8z2fnydpWdIdki4e9iVJO9P5lldWVqZLsZmVq8wm+4O6LPnhD48PJlBaNyZNVGVA0YBxA7NLkt4JLAK/2zN6Lo2Yvwb8nqSXDfpuROyJiMWIWJydnZ02zdZmLtqonzLv+h4nSPketIGqDCiHgdN7Pp8GHOmfSdJbgF3ARRHxzOr4iDiSvj4M3A68psjEWsu5e516KrPJ/jhByvegDVRlQLkTOEvSGZI2ApcAx7XWSutNPkoSTJ7sGX+KpJPT91uANwAPlJZya58Ke2i1Ecpssj8oeG3cCCeddPw434M2VGV3ykfEs5KuAG4B1gPXR8T9kq4meZjLfpIirhcAn5EEcCht0fVK4KOSjpIExWsiwgHFJufudeqrrLu+V5exa1fyv8/NHQsc/eM6VCE/DncOaQbuANRshCa08jKrD3evYzY1BxQzcPc6ZjlwQDFb5e51LE8dbIbu7uvNzPLWf9f9ajN0aPWFinMokxj3yqODVypmndbRZujOoYxr3CuPjl6pmHVaR5uhu9nwuMZtXurmqGbd07Lj3s2GizLulUdHr1TMOq2jzdAdUMY1bmd1ZXZuZ2b10NFm6A4o4xr3yqNJVypuPGCWnw42Q3dAGde4Vx5NuVJxb7tmxenIxZor5S3RskpEs9zt2zdZJ5GDngTZsKc+Zq2Ud0CxxLp1Sc6kn5Rk2c26bJqg0IKLNbfysvG48YDZcNPcqNihlp4OKJZoUuMBs7JNExQ6dLFWaUCRdIGkByUdlHTlgOknS7opnf5VSQs9065Kxz8o6a1lpruVmtJ4oEk6UhHbCdMEhS5drEVEJQPJUxofAs4ENgLfAM7um+dy4Lr0/SXATen7s9P5TwbOSH9n/VrL3LZtW5iVYu/eiJmZiKRmKhlmZpLx1jzT/p9790bMz0dIyWvD9gOSp+iueV6vModyLnAwIh6OiB8CnwK2982zHbgxfX8zcJ6SZwFvBz4VEc9ExCPAwfT3zOqhqZ0DOlc12LQ5+DLvSanwP8wUUCTdJunCvnF7plz2qcDjPZ8Pp+MGzhMRzwLfAzZn/K5ZdZpYEet7kUZrwo2KFf+HWXMoZwAfkPTBnnFrNiFbgwaM62+3OmyeLN9NfkDaKWlZ0vLKysqYSTSbUBMrYpuaq7JjKv4PswaU7wLnAS+V9OeSXpTDsg8Dp/d8Pg04MmweSRuAFwFPZfwuABGxJyIWI2JxdnY2h2SbZdDEitgm5qrseBX/h1kDiiLi2Yi4HPgz4L8DL5ly2XcCZ0k6Q9JGkkr3/X3z7Ad2pO/fDnwxrSDaD1yStgI7AzgL+NqU6THLTxNbzTUxV2XHq/g/zBpQrlt9ExE3AJcCfznNgtM6kSuAW4D/BXw6Iu6XdLWki9LZPg5slnQQeB9wZfrd+4FPAw8A/wX4pxHx3DTpMctdE8rcezUxV2XHq/g/dNcrZnbMpP1VWX0U8B+6L68BHFDMzMbnvrys27pyP0VX1tMaYUPVCTDLXX/PsKtt8aFdxTddWU9rDBd5Wfu0oLvwTLqynlY5F3lZd3XlfoqurGdduHhxTQ4o1j5duZ+iK+tZB+6WJhMHFGufrtxP0ZX1rAN3S5OJA4q1TxPvUp9EV9azDly8mIkr5c3M1tLxBhCulDczy4uLFzNxQDEzW4uLFzPxjY1mZlksLTmArME5FDMzy4UDirWDbzozq5yLvKz53KeVWS04h2LN55vOzGqhkoAiaZOkWyUdSF9PGTDPOZL+StL9ku6R9I96pt0g6RFJd6fDOeWugdWKbzozq4WqcihXArdFxFnAbennfk8D746IVwEXAL8n6cU9098fEeekw93FJ9lqy31amdVCVQFlO3Bj+v5G4OL+GSLimxFxIH1/BHgSmC0thdYcvunMrBaqCigvjYgnANLXl4yaWdK5wEbgoZ7Ru9OisGslnVxcUq32fNOZWS0U1peXpC8APzlg0i7gxoh4cc+834mIE+pR0mlbgduBHRFxR8+4/0sSZPYAD0XE1UO+vxPYCTA3N7ftsUH98ZiZ2VCV9+UVEW+JiFcPGD4LfCsNCqvB4clBvyHpJ4D/DPyr1WCS/vYTkXgG+BPg3BHp2BMRixGxODvrEjOzTvJ9SqWoqshrP7Ajfb8D+Gz/DJI2Av8R+EREfKZv2mowEkn9y32FptbMmssPxypNVQHlGuB8SQeA89PPSFqU9LF0nl8FfgG4dEDz4H2S7gXuBbYAHy43+WZWmXFzG75PqTR+HoqZNUd/rwiQtOgb1Qhj3bokZ9JPgqNHi0lny1Reh2JmlrtJchu+T6k0Dihm1hyT9Irg+5RK44BiZs0xSW7D9ymVxgHFzJpj0tzG0lLy7PejR5NXB5NCOKCYWXM4t1FrDijWXL5ZrZuc26gtBxRrJt+s1ky+CGg1BxQ7pkkHu29Wa566XgQ0ab+vOd/YaIlJbhirkm9Wa56FhSSI9JufT4quqtC0/b4iWW9sdECxRB0P9lGall6r50WA96NMfKe8jadpj9H1zWrNU8c71pu239ecA4ol6niwj+Lmo81Tx4uApu33NeeAYok6HuxrcfPRZqnjRUAT9/sac0ApW11blNTxYLf2qdtFgPf7XLlSvkxuUWJmDeRK+TryvRNm1mKVBBRJmyTdKulA+nrKkPme63la4/6e8WdI+mr6/ZvSxwXXn1uUmFmLVZVDuRK4LSLOAm5LPw/yg4g4Jx0u6hn/b4Br0+9/B7is2OTmxC1KzKzFqgoo24Eb0/c3Ahdn/aIkAW8Gbp7k+5VyixIza7GqAspLI+IJgPT1JUPme56kZUl3SFoNGpuB70bEs+nnw8CpwxYkaWf6G8srKyt5pX8yblFiZnmpYYvRDUX9sKQvAD85YNI4NdBzEXFE0pnAFyXdC/y/AfMNbaoWEXuAPZC08hpj2cVYWnIAMbPp9LcYXe1oEyo9vxQWUCLiLcOmSfqWpK0R8YSkrcCTQ37jSPr6sKTbgdcAfwa8WNKGNJdyGnAk9xUwM6urUS1GKwwoVRV57Qd2pO93AJ/tn0HSKZJOTt9vAd4APBDJjTNfAt4+6vtmZq1V0xajVQWUa4DzJR0Azk8/I2lR0sfSeV4JLEv6BkkAuSYiHkinfQB4n6SDJHUqHy819WZmVappi9HCirxGiYhvA+cNGL8M/Eb6/ivAzwz5/sPAuUWm0cystnbvHtzrRsUtRn2nvLVPDVu/mOWqpi1GK8mhmBWmpq1fzHJXwxajzqFYu7i/NLPKOKBYu9S09YtZoWpSzOuAYu1Sh9YvNTm4rSNWi3kfewwijhXzVrDfOaBYu1TdX1qNDm7riBoV8zqgWLtU3fqlRge3dUSNinndysvap8rWLzU6uK0j5uaSnPCg8SVzDsUsT3Wow7FuqbqYt4cDilmeanRwW0dUXczbw0VeZnlaPYh37UqKuebmkmBSsxvQrGVqcpOjcyhFcLPRbltagkcfhaNHk9caHOg2go/X3DiHkjd3/WHWHD5ec6Xk8SLdsLi4GMvLy8UuZGFhcIuL+fnkatXM6sPHayaS7oqIxbXmc5HXuNbKHrvZqFlz+HjNlQPKOLLcBe1mo2bN4eM1V5UEFEmbJN0q6UD6esqAed4k6e6e4W8lXZxOu0HSIz3Tzikl4VnugnazUbPm8PGaq6pyKFcCt0XEWcBt6efjRMSXIuKciDgHeDPwNPCXPbO8f3V6RNxdSqqzZI9r1CbcrFRNbC3l4zVXlVTKS3oQeGNEPCFpK3B7RLxixPw7gV+MiKX08w3AX0TEzeMsd+pKeVfgmQ3W31oKkit9n5xboe6V8i+NiCcA0teXrDH/JcAn+8btlnSPpGslnTzsi5J2SlqWtLyysjJdqp09NhvMnWIaBQYUSV+QdN+AYfuYv7MV+Bnglp7RVwF/F/g5YBPwgWHfj4g9EbEYEYuzs7Pjr0hvNn7XLtixw9njNqmqmKaJxUOjuLWUAURE6QPwILA1fb8VeHDEvP8C2DNi+htJir/WXO62bdtiLHv3RszMRCRtupJhZiYZb81X1f/bxv1qfv749Vkd5uerTpnlAFiODOfYqoq89gM70vc7gM+OmPcd9BV3pbkWJAm4GLivgDQ6G992Vf2/bdyvXBxsVFeHcg1wvqQDwPnpZyQtSvrY6kySFoDTgf/a9/19ku4F7gW2AB8uJJXOxrdbVf9vG/erNreWalvxZIHc9coobtXVblX9v96vmsOt14D6t/JqBmfj262q/9f7VXO0sXiyQA4oozQ9G++s+mhV/b9N36/arP+YGZSThGYXTxbIRV5t5ay62XgGHTNS0l5tkPn5zjw8zUVeXVdkVt05H2ujQcdMRBJUBhnUOWzHOaC0VVEtibL0uGzWRMOOjYgkNzKI61OO44DSVkV1y93USkrnqmwtw44NKSnaGpZTcX3KjzmgtFVRLYmaeA+Fc1WW5YJiWNCISC6Y/OyUNTmgtFVRLYmaeFA1NVdl+ch6QbG0NLwC/tAhN/fOwAGlzZaWkhvljh5NXvNojdLEg6qJuSrLzzgXFMPqSubm3Nw7AwcUG08TD6om5qosP+NcUKx1wVTERVqLOKDY+Jp2UDUxV2X5GeeCookXTDXigGLt55NEt417QdG0C6Ya2VB1AsxKsbTkE0NXrf7vu3YlxVxzc525w71szqFMyvc1JLwdbJg67RvOdZTCAWUSvq8h0dTtUKcTXdMN25ZN3TdsOlke65j3APwKcD9wFFgcMd8FJI8LPghc2TP+DOCrwAHgJmBjluWO/QjgYdr2uNO9e5O0S8lr1kfRVr0dJkl3Gx+/W5VR2zLPfWPS/dNyQ8ZHAFcVUF4JvAK4fVhAAdYDDwFnAhuBbwBnp9M+DVySvr8OeG+W5eYWUKTBB4uUz++XaZoTbJXbYdJ0Vx0E22TYtly/fvD4SfaNPC8A2hyYCl63WgeUHy98dEB5PXBLz+er0kHAXwMbBs03anAOZYBp1qXK7TDpstt0MVC1Ydty1DDuvpHXPtbmnGkJ65Y1oNS5DuVU4PGez4fTcZuB70bEs33jy9Om+xqmuYu8yu0wabp9k2N+xt1mk+wbefVy0Obud2q0boUFFElfkHTfgGF71p8YMC5GjB+Wjp2SliUtr6ysZFz0Gtp0X8M0J9gqt8Ok6W7TxUDVBm3LYSbdN/K6AGhz9zt1Wrcs2ZiiBppa5NUmTS0KmCbdbS5LL1vvthxWdzJNEWhe+2ebiqn7lbButKAOZQPwMEmLrtVK+Vel0z7D8ZXyl2dZngPKEE09wTY13W1V1MVJHv9zUy+csqhRHUpVgeSXSeo+ngG+tZrDAH4K+FzPfBcC3yRp7bWrZ/yZwNdImhN/Bjg5y3IdUMwKVucgX+e0TasmrbyUzNsNi4uLsby8XHUyzMwaRdJdEbG41nx1buVlZmYN4oBiZma5cEAxM7NcOKCYmVkuHFDMzCwXnWrlJWkFeGzMr20huZGyS7q4ztDN9fY6d8O06zwfEbNrzdSpgDIJSctZmsu1SRfXGbq53l7nbihrnV3kZWZmuXBAMTOzXDigrG1P1QmoQBfXGbq53l7nbihlnV2HYmZmuXAOxczMcuGAkpJ0gaQHJR2UdOWA6SdLuimd/lVJC+WnMl8Z1vl9kh6QdI+k2yTNV5HOPK21zj3zvV1SSGpFa6As6y3pV9P/+35Jf1p2GvOWYf+ek/QlSV9P9/ELq0hnniRdL+lJSfcNmS5Jf5Buk3skvTbXBGTpkrjtA7CepIv8Mzn27JWz++a5HLgufX8JcFPV6S5hnd8EzKTv39uFdU7neyHwZeAOhjyvp0lDxv/6LODrwCnp55dUne4S1nkP8N70/dnAo1WnO4f1/gXgtcB9Q6ZfCHye5EGFrwO+mufynUNJnAscjIiHI+KHwKeA/kcVbwduTN/fDJwnadDjiJtizXWOiC9FxOrDqu8ATis5jXnL8j8D/BbwO8Dflpm4AmVZ738CfCQivgMQEU+WnMa8ZVnnAH4iff8i4EiJ6StERHwZeGrELNuBT0TiDuDFkrbmtXwHlMSpwOM9nw+n4wbOExHPAt8DNpeSumJkWedel5Fc2TTZmuss6TXA6RHxF2UmrGBZ/uuXAy+X9D8k3SHpgtJSV4ws6/wh4J2SDgOfA/5ZOUmr1LjH/Vg25PVDDTcop9Hf/C3LPE2SeX0kvRNYBH6x0BQVb+Q6S1oHXAtcWlaCSpLlv95AUuz1RpKc6H+T9OqI+G7BaStKlnV+B3BDRPxbSa8H/n26zkeLT15lCj2POYeSOAyc3vP5NE7M/v54HkkbSLLIo7KWdZdlnZH0FmAXcFFEPFNS2oqy1jq/EHg1cLukR0nKmPe3oGI+6/792Yj4UUQ8AjxIEmCaKss6XwZ8GiAi/gp4HkmfV22W6biflANK4k7gLElnSNpIUum+v2+e/cCO9P3bgS9GWsvVUGuuc1r881GSYNL0MnVYY50j4nsRsSUiFiJigaTe6KKIaPpzo7Ps3/+JpBEGkraQFIE9XGoq85VlnQ8B5wFIeiVJQFkpNZXl2w+8O23t9TrgexHxRF4/7iIvkjoRSVcAt5C0Drk+Iu6XdDWwHBH7gY+TZIkPkuRMLqkuxdPLuM6/C7wA+Eza/uBQRFxUWaKnlHGdWyfjet8C/ANJDwDPAe+PiG9Xl+rpZFznfwn8saTfJCn2ubThF4lI+iRJseWWtG7og8BJABFxHUld0YXAQeBp4NdmB95gAAAA8klEQVRzXX7Dt5+ZmdWEi7zMzCwXDihmZpYLBxQzM8uFA4qZmeXCAcXMzHLhgGJmZrlwQDEzs1w4oJhVSNLPpc+leJ6kv5M+i+TVVafLbBK+sdGsYpI+TNLtx/OBwxHx2xUnyWwiDihmFUv7mrqT5PkrPx8Rz1WcJLOJuMjLrHqbSPpMeyFJTsWskZxDMauYpP0kTxQ8A9gaEVdUnCSzibi3YbMKSXo38GxE/Kmk9cBXJL05Ir5YddrMxuUcipmZ5cJ1KGZmlgsHFDMzy4UDipmZ5cIBxczMcuGAYmZmuXBAMTOzXDigmJlZLhxQzMwsF/8fuQAK78VbanMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1590a5605c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.plot(x,z,'ro')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们从导入`ambhas.krige`库开始。我们使用其名称为`OK`的函数Oridinary Kriging来创建一个Kriging类的实例。为了查看变差函数，我们使用`variogram`方法。有两种类型的变量函数可用，第一个是原始变差图，即数据的平均不可执行，第二个是在均匀的滞后距离区间上的平均。让我们先看看原始的变差图是怎么样的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amhas import krige\n",
    "foo = krige.OK(x,y,z)\n",
    "D,G = foo.variogram(var_type='scattered')\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(D,G,'.')\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "原始变差图如图10.4所示。这个变差图非常离散，并不能提供对数据行为的任何了解。为了观察数据的行为，我们使用平均变差函数。\n",
    "\n",
    "平均变差函数的估计(或简单的称为变差函数)是通过将`averaged`作为`variogram`方法的参数`var_type`的输入来实现。变差图如图10.5所示。该变差图看上去趋于水平。这意味着行为在任何滞后距离都是相同的，或变量(z)在任何尺度上都有相同的变化。这可能是因为变量就是如此，或者我们的抽样太差，以至于我们无法捕捉变量的行为。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DE,GE = foo.variogram(var_type='averaged',min_lag=0.025)\n",
    "\n",
    "plt.clf()\n",
    "plot.plot(DE,GE,'ro')\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.ylim((0,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "让我们尝试另一个变量，它在比我们的采样更大规模上有变异性。我们保留相同的采样，并以下面的方式改变信号的变异性：\n",
    "<center>$z(x)=cos(2x/0.25)\\quad\\quad\\quad\\quad\\quad\\quad(10.2)$</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X+QXeV93/H3Z3cRtoDYYiXbBNCKpEpq4nawuyFOM9PGBtsyfyDakgS8AiHjKmZLy8Stx2TUGTo4mtrNZLDSsWRkWyDDxoBJU5TGLuGX67YxlGWCMYLByLIWVFEjJPyDyEFI+vaPc250dXXu3bu7997z6/OaubP3POe5d58j3bvf8/xWRGBmZrZQQ3kXwMzMqsEBxczMesIBxczMesIBxczMesIBxczMesIBxczMesIBxczMesIBxczMesIBxczMemIk7wIM0tKlS2PFihV5F8PMrFSeeOKJVyJi2Wz5ahVQVqxYwfT0dN7FMDMrFUkz3eRzk5eZmfWEA4qZmfWEA4qZmfWEA4qZmfWEA4qZmfVErgFF0jZJL0t6us15SfpjSbskPSXpPU3n1kp6Pn2sHVypzazwJidhaAik5HHGGTA1lXepKi/vGsrtwKoO5z8MrEwf64EtAJLOBG4Cfg24ELhJ0pK+ltTMim1qClasSALIli3QvBvta6/BNdc4qPRZrgElIr4FHOyQZTXwlUg8CrxV0lnAh4AHIuJgRLwKPEDnwDR/jQ/p0FDy0x9Is2KZmoKlS2HNGpjpMF3iyBHYsGFw5aqhvGsoszkbeLHpeG+a1i79JJLWS5qWNL1///65/fapKVi/PvmQRiQ/16xJ7oBGRpJqtZnlp/EdPXCgu/wvvNDf8tRc0QOKMtKiQ/rJiRFbI2I8IsaXLZt15YATbdgAhw5lnzt6NKlWX3zx3N7TzHqn03c0y/Ll/SuLFT6g7AXObTo+B9jXIb23urmbeeghN4OZDUprE3SnJq5WIyOwcWPnPJOTST63QsxL0QPKDuDqdLTXe4EfR8RLwP3AByUtSTvjP5im9Va3dzNulzXrv6wmaGU1VmQ4/XS4/XaYmGifZ3IyaXU4ejQ5brRCOKh0Le9hw18Fvg38sqS9kq6V9HFJH0+zfB3YDewCvghMAkTEQeDTwOPp4+Y0rbc2boTFi2fPNzPjTnuzfstq3orIDiqjo3Dnncn5CPjpTzsHE4CtW7PTt2zx97pLisjseqik8fHxmPNqw1NTyQe526r14sXJB3O2D6+Zzc3Q0IlDgZuNjSVN1MuXJzeC8/n+dart1Px7LemJiBifNZ8DyhxcfHHSZzKb4WHYvr22Hz6zvmjXZzI2Bnv2LPz9R0aON3dl6dXvKaFuA0rR+1CK5cEHk2r02Fjnu5mjR5O2XleTzXonqwl68eLZO9q7tX595/MecjwrB5S5mphI7lKOHUsCSzuHDsHatQ4qZr0yMZE0OzVu6MbGetsMtXkzXHdd+/MecjwrB5SFmK3T3jUVs95qvqHbs6f3zcqbNyetEP2sCVWYA8pCNO6Yhofb5zl0yMOKzbpRlGWOOtWEilLGgnKnfC80xse3m7ErJXdUZpYt6ztUtJFVZShjn3iUV4a+BRRIPmxr12aPEqnx6BCzrpx+OvzN35ycXqTvTrtRZqOj8MorAy/OIHmU16BNTCRDhd32ajY3k5PZwQSKNbKqXVkOHHDTV8oBpZf6PQrFrGqmppKZ6O0UaWRVp7K4nxRwQOm9dqNQ3JlndqKpKbj66s55ilS771SWItWkcuQ+lEGYmoJ16+CNN46nnXIK3Habay9WX+36TRqGhjrPXM/D0qXZe68Uqa+nD9yHUiQ33HBiMIHk+IYb8imPWd6mpjoHE4Df/d3BlGUuNm1yP2kHDiiD0G43uW53mTOrmtn6HE47LZlkWDTuJ+3IAcXMBm+2Podbbx1MOeYjq5/UG3MBDiiDMTqanT405M55q6dOI6auu65cd/zemOvv5L3B1ipJz0naJenGjPO3SHoyfXxP0o+azh1tOrdjsCWfo02bYNGik9OPHfNaX1ZPWevgSUkwKWJTVyftNuZql15huQUUScPA54EPA+cDV0o6vzlPRPxeRFwQERcA/xn4L02nf9Y4FxGXDqzg8zExAdu2Za/55bW+rI6y+iLuuKN8wQTaj0Q7erR2N4t51lAuBHZFxO6IOAzcBazukP9K4KsDKVk/TEy0X89rZsZzU6z6WudiQX9XDh6UTovDrllTq6avPAPK2cCLTcd707STSBoDzgMebkp+k6RpSY9Kuqx/xeyhTu3GMzPw0Y86qFg1NRZWnJlJtvGdmalOc+9sG3PVaE/6PANK1paH7WZZXgHcGxHNdcvl6USbjwCfk/SLmb9EWp8Gnun9+/cvrMQLNdv+KYcPe26KVdOGDSevxl2V5t7ZNuaCalxnF/IMKHuBc5uOzwH2tcl7BS3NXRGxL/25G/gm8O6sF0bE1ogYj4jxZcuWLbTMC9PcbtyO56ZYFbUbJlyVJUtm6/upynXOIs+A8jiwUtJ5khaRBI2TRmtJ+mVgCfDtprQlkk5Nny8FfgN4ZiClXqjGGHazOmnX3FukxR8XaqjDn9OaTBHILaBExBHgeuB+4FngnojYKelmSc2jtq4E7ooTFx17JzAt6TvAI8BnIqIcAaWh3dyUdulmZZbV3Fu1JUs6LRVTk+3AvThkXrIWjGwYG0u+aGUd9WKWZWoq6Ut44YWkZlLFz/jkZDLLv92IzpIuIukdGzMUKqDA8S/YzEwyFr/5/6ImW4uaVdLQ0Inf54aSbgfu1YbLoNGfMjZ28oevKiNgzOqoDn1GGRxQiqDqI2DM6qbdFIHXXqt0P4oDShHU9G7GrLIaUwRaB9kcOFDpznkHlCKowwgYs7qZmEh2pWxV4eZsB5Qi8KY9ZtVUs+bskbwLYKmJCQcQs6pZvjwZxZmVXkGuoZiZ9UvNmrMdUMzM+qVmzdkOKEXWun9ERUeGmFVa1h70FeWAUlRZ+0dcdVWtNusxs3JxQCmqrP0jIuALX3BNxYrFNWlLOaAUVbthhRGVHcNuJVTlnRhtzhxQiqrTsMKKjmG3EqryToz9VsGanQNKUW3cmIwKyXLmmYMti1k7NZu41zMVrdk5oBTVxAR8/OPZ537609J/8KwivA7d/FS0ZpdrQJG0StJzknZJujHj/DWS9kt6Mn18rOncWknPp4+1gy35gGzenL2D4+HDcMMNgy+PWauaTdzrmYrW7HILKJKGgc8DHwbOB66UdH5G1rsj4oL08aX0tWcCNwG/BlwI3CRpyYCKPlgHD2anHzjgWorlr2YT93qmojW7PGsoFwK7ImJ3RBwG7gJWd/naDwEPRMTBiHgVeABY1ady5qvTB6zk1WOriBpN3OuZitbs8gwoZwMvNh3vTdNa/QtJT0m6V9K5c3xt+XX6gJW8emxWW601u9FRePObk8nLJR7xlWdAyRrC1LoJ858DKyLiHwIPAtvn8Noko7Re0rSk6f3798+7sLmZmMjuR4HSV4/Naq1Rs7vjDvjZz5Jm7JKP+MozoOwFzm06PgfY15whIg5ExOvp4ReBf9Tta5veY2tEjEfE+LJly3pS8IHbtKmS1WMzo1IjvvIMKI8DKyWdJ2kRcAWwozmDpLOaDi8Fnk2f3w98UNKStDP+g2laNbnj06y6KjTiK7cNtiLiiKTrSQLBMLAtInZKuhmYjogdwL+RdClwBDgIXJO+9qCkT5MEJYCbI6LNcKiK8AZcZtVUoU24FJHZ9VBJ4+PjMT09nXcxzMyOa8yab272Wry4UK0Qkp6IiPHZ8nmmfFlNTsLISNIENjLiZe3NyqpCTdreU76MJidhy5bjx0ePHj/evDmfMpnZ/FWkSds1lDLaunVu6WZmA+CAUkZHj84t3cxsABxQymh4uP25Ek6GsoKr4L4dhVfSf3MHlDJav779Oe87b700NQXr1p24b8e6daX5A1dKJd4rxcOGy6q1Y77VnXdWopPPcrZ0abIkSKvRUXjllcGXpw5WrMielzI2lizVkoNuhw07oJTZ0FByB5PFX3jrhXa7hkL7z54tTLvvtZSs6JwDz0Opg04zabPuKs3mogRNLJVU4r1SHFDKzItDWr80+k7aabcCti1cifdKcUAps4kJOO207HP+wttCbNgAb7zR/vymTYMrS92UeK8UB5Syu/VWOOWUE9NOOcVfeFuY2Va69YCP/irpXikOKGU3MQG33XbiOkC33eYvvC1Mp/b6sbHBlaPuSrZXigNKFXhPb+u1Sy7JTl+0qBRt+ZVRsr1SHFDM7ERTU7B9+8npp50G27b5hmWQSjbiywGlikq6bIMVRFYzCySTHB1MBqtkI75yDSiSVkl6TtIuSTdmnP+EpGckPSXpIUljTeeOSnoyfexofW1tlXjZBiuIkjWzVFrJ9krJbaa8pGHge8AHgL0k2/leGRHPNOV5H/BYRBySdB3wmxHxO+m51yLi9Ln8zsrNlM9SwGUbrGT8GbIWZZgpfyGwKyJ2R8Rh4C5gdXOGiHgkIhp170eBcwZcxvLx3aUtVMmaWaw48gwoZwMvNh3vTdPauRb4RtPxmyRNS3pU0mXtXiRpfZpvev/+/QsrcRmUrBPPCqTR93bVVclEutHRUjSzWHHkGVCyVp3LbH+TtAYYB/6wKXl5WgX7CPA5Sb+Y9dqI2BoR4xExvmzZsoWWufh8d2nz0dr3duBAMqHujjs8FN26lmdA2Quc23R8DrCvNZOki4ENwKUR8XojPSL2pT93A98E3t3PwpZGyTrxrCBKNoGu1go8ijPPTvkRkk75i4D/S9Ip/5GI2NmU593AvcCqiHi+KX0JcCgiXpe0FPg2sLq5Qz9LLTrlzeajgEumW4ZGTbI5+C9e3PebxsJ3ykfEEeB64H7gWeCeiNgp6WZJl6bZ/hA4Hfhay/DgdwLTkr4DPAJ8ZrZgUlsFvpuxAnHfWzkUvCbpDbaqLKe7GSshf1bKIaeaZOFrKDYABb+bsQJx31s5FLwm6YBSZZ6TYnPhRUaLr+CjOB1QqqzgdzNmNkcFr0k6oFRZwe9mzGweClyTdECpsoLfzZhZtYzkXQDrs4kJBxAzGwjXUMzMrCccUMzMrCccUMzMrCccUMzMrCccUMzMrCccUMzMrCccUOrIKxCbWR94HkrdtK4qOzOTHIPnq5jZgriGUjdegdjM+iTXgCJplaTnJO2SdGPG+VMl3Z2ef0zSiqZzv5+mPyfpQ4Msd6l5BWIz65PcAoqkYeDzwIeB84ErJZ3fku1a4NWI+HvALcBn09eeD1wB/AqwCticvp/Npt1Kw0ND7ksxswXJs4ZyIbArInZHxGHgLmB1S57VwPb0+b3ARZKUpt8VEa9HxA+AXen72WyyViAGOHo06UtxUKmmyUkYGUkWCR0ZSY7NeqyrgCLpIUmXtKRtXeDvPht4sel4b5qWmSfdg/7HwGiXr7UsjRWIhzMqdO5LqabJSdiyJblpgOTnli0OKtZz3dZQzgM+JemmprRZ9xeehTLSWjdLbpenm9cmbyCtlzQtaXr//v1zLGJFTUy033/afSnVc+ut2elbF3pPaHaibgPKj4CLgLdL+nNJb+nB794LnNt0fA6wr10eSSPAW4CDXb4WgIjYGhHjETG+bNmyHhS7IrybYz1MTbW/eWjUWKyacphv1m1AUUQciYhJ4E+B/wW8bYG/+3FgpaTzJC0i6WTf0ZJnB7A2fX458HBERJp+RToK7DxgJfB/FlieevFujvXQqQkzq9nTqqEx32xmBiKOzzfrc1DpNqB8ofEkIm4HrgH+ciG/OO0TuR64H3gWuCcidkq6WdKlabYvA6OSdgGfAG5MX7sTuAd4BvjvwL+KCN9uzYV3c6yHTk2YjQmtVj05zTdTcsNfD+Pj4zE9PZ13McwGZ8WK5O601WmnwWuvDbw4NiBDQ0nNpJXUvgm0A0lPRMSs/eaeKW9WVVNT2UFj8eL2HfVWDe36Qs88s6+/1gHFrIoabegHDpyYPjrqps062LgRFi06Of0nP+lrP4oDih3nVYirI6sNHeD00x1M6mBiAs444+T0N97oaz+KVxu2hFchrhav2WYHD2an9/Ez4BqKJbwKcbV4npHl8BlwQLGE72iro1NnvOcZ1UcOc80cUCzhO9pqmJqCdevcGW+5zDXzPBRLtPahQHI34z9C5bJ06cnBBJKA8sorgy+PVYLnodjceOZ8NWQFk07pZj3kUV523MSEA4iZzZtrKGZVMjo6t3SzHnJAMauSTZtOniG9aFGSbtZnDihmVTIxAdu2ndgXtm2bmzJtINyHYlY17guznLiGYlZ2XoPNCsI1FLMy8xpsViC51FAknSnpAUnPpz+XZOS5QNK3Je2U9JSk32k6d7ukH0h6Mn1cMNgrqBnfAReX12CzAsmryetG4KGIWAk8lB63OgRcHRG/AqwCPifprU3nPxkRF6SPJ/tf5JrKaW9q65LXYLMCySugrAa2p8+3A5e1ZoiI70XE8+nzfcDLwLKBldASvgMuNq/BZgWSV0B5e0S8BJD+fFunzJIuBBYB329K3pg2hd0i6dT+FbXmfAdcbDmsKGvWTt8CiqQHJT2d8Vg9x/c5C7gDWBcRx9Lk3wf+PvCrwJnApzq8fr2kaUnT+/fvn+fV1JjvgIvNa7BZgfQtoETExRHxrozHfcAP00DRCBgvZ72HpJ8D/gL49xHxaNN7vxSJ14HbgAs7lGNrRIxHxPiyZW4xmzPfARdXY7DEVVclx3fcAXv2OJhYbvJq8toBrE2frwXua80gaRHwZ8BXIuJrLecawUgk/S9P97W0deY74GLyYAkroFz2Q5E0CtwDLAdeAH4rIg5KGgc+HhEfk7SGpPaxs+ml10TEk5IeJumgF/Bk+pqMLepO5P1QrDJWrEiCSKuxsaSWYtZD3e6H4g22zMpoaCipmbSS4Nixk9PNFsAbbJlVmQdLWAE5oJiVkQdLWAE5oJiVkQdLWAF5cUizsvIy9VYwrqHY/HjBSDNr4RqKzZ2XTDezDK6h2Nx5wUgzy+CAYnPnBSPNLIMDis2d50CYWQYHFJs7z4EwswwOKDZ3ngNhZhk8ysvmx3MgzKyFayjWG56XYlZ7Dii2cN6boz8cpK1kHFBs4TwvpfccpK2EHFBs4TwvpfccpK2Ecgkoks6U9ICk59OfS9rkOyrpyfSxoyn9PEmPpa+/O90u2PLieSm95yBtJZRXDeVG4KGIWAk8lB5n+VlEXJA+Lm1K/yxwS/r6V4Fr+1tc68jzUnrPQdpKKK+AshrYnj7fDlzW7QslCXg/cO98Xm994HkpvecgbSWUV0B5e0S8BJD+fFubfG+SNC3pUUmNoDEK/CgijqTHe4Gz+1tcm9XEBOzZk+xnvmePg8lCOUhbCfVtYqOkB4F3ZJyaS6/i8ojYJ+kXgIclfRf4SUa+6FCO9cB6gOVuLrCim5pKOt5feCFp3tq40UHESqNvASUiLm53TtIPJZ0VES9JOgt4uc177Et/7pb0TeDdwJ8Cb5U0ktZSzgH2dSjHVmArwPj4eNvAY5Y77zNjJZdXk9cOYG36fC1wX2sGSUsknZo+Xwr8BvBMRATwCHB5p9eblY6HClvJ5RVQPgN8QNLzwAfSYySNS/pSmuedwLSk75AEkM9ExDPpuU8Bn5C0i6RP5csDLb1ZP3iosJVcLotDRsQB4KKM9GngY+nzvwL+QZvX7wYu7GcZzQZu+fKkmSsr3awEPFPerCg8VNhKzgHFrCg8VNhKzgHF+ssr5s6N5/NYiXmDLesfD4M1qxXXUKx/PAzWrFYcUKx/PAzWrFYcUKx/vGKuWa04oFj/eBisWa04oFj/eBisWa04oFh/tQ6DBQ8jNqsoDxu2wfEwYrNKcw3FBsfDiM0qzQHFBsfDiM0qzQHFBsfDiM0qzQHFBsfDiM0qzQHFBsfDiM0qLZeAIulMSQ9Iej79uSQjz/skPdn0+FtJl6Xnbpf0g6ZzFwz+KmxevJquWWXlVUO5EXgoIlYCD6XHJ4iIRyLigoi4AHg/cAj4y6Ysn2ycj4gnB1Jqs/nwEv5WE3kFlNXA9vT5duCyWfJfDnwjIg7Nks/Kqqp/dBtzb2ZmIOL43JuqXJ9Zk7wCytsj4iWA9OfbZsl/BfDVlrSNkp6SdIukU9u9UNJ6SdOSpvfv37+wUlt/VPmPrufeWI0oIvrzxtKDwDsyTm0AtkfEW5vyvhoRJ/WjpOfOAp4Cfj4i3mhK+3/AImAr8P2IuHm2Mo2Pj8f09PScr8X6bMWKJIi0Ghs7vlxLWQ0NJUGylZT0I5mVgKQnImJ8tnx9W3olIi5ud07SDyWdFREvpcHh5Q5v9dvAnzWCSfreL6VPX5d0G/DvelJoy0e7iY1ZQaZsli/Pvg7PvbEKyqvJawewNn2+FrivQ94raWnuSoMQkkTS//J0H8pog9Lpj+vk5ODK0SuTkzAyktRCXnwxqaU089wbq6i8AspngA9Ieh74QHqMpHFJX2pkkrQCOBf4Hy2vn5L0XeC7wFLgDwZQZuuXTn9ct24dXDl6YXIStmyBo0eT42PHksdpp3nujVVe3/pQish9KAUmtT9Xps/oyMjxYNJseBiOHBl8ecx6oNs+FM+Ut2IYHm5/bmSkPE1fWcGkU7pZhTigWDE09kXJcvRo0oxU5KDSmEfTTqeAaVYRDihWDJs3w3XXdf7DW9T+lOZ5NO10CphmFeGAYsWxeXPnfoaiNhtlTV5sGB5OAuXmzYMtk1kOvAWwFc/wcPuO7SJqN49Gcke81YprKFY87ZqHitZs1Og3aTcKzZMXrWZcQ7HiaTQPbd2a1FSGh5NgUqRmo0a/SbumLk9etBpyDcWKqdGfEpH8LFowufrq9sHEkxetplxDMZuLqSn46EfbL+wolX9BS7N5cg3FbC42bIDDh9ufd7+J1ZgDipXT1BQsXZrUCKTk+SD2T2k3oqvB/SZWYw4oVj6NZqcDB46nHTgA69b1P6h0qoGMjrrfxGrNAcXKp12z0xtv9H8nxI0bYdGik9NPOQU2berv7zYrOAcUK59OzU4zM/3dk35iArZtS2ojDaOjcNttrp1Y7TmgWPnM1vE9MwNr1vSvX2ViAl55JRnSHJE8dzAxyyegSPotSTslHZPUdo19SaskPSdpl6Qbm9LPk/SYpOcl3S0pow3CKqtds1OrAweSyYeD6Kw3s9xqKE8D/xz4VrsMkoaBzwMfBs4HrpR0fnr6s8AtEbESeBW4tr/FtULJanZq59ChE/tVGsulDA31t2nMrIZyCSgR8WxEPDdLtguBXRGxOyIOA3cBq9N95N8P3Jvm206yr7zVSXOz09hY57yNPpfmZeYjkp+uwZj1TJH7UM4GXmw63pumjQI/iogjLelWVxs3JmtntdPoc8laZr61BmNm89a3pVckPQi8I+PUhoi4r5u3yEiLDuntyrEeWA+w3LOYq6nRIX7DDSfOTYETF2lsNzpstsmKZtaVvtVQIuLiiHhXxqObYAJJzePcpuNzgH3AK8BbJY20pLcrx9aIGI+I8WXLls3nUqwMGk1gd96ZNIFJJy/S2O6GwjcaZj1R5Cavx4GV6YiuRcAVwI6ICOAR4PI031qg2yBlVTcxkSzOeOxY8rN5OG9W05iXmTfrmbyGDf8zSXuBXwf+QtL9afrPS/o6QNpHcj1wP/AscE9E7Ezf4lPAJyTtIulT+fKgr8FKaGIiqbG0q8GY2YIo2u02V0Hj4+MxPT2ddzHMzEpF0hMR0XbOYEORm7zMzKxEHFDMzKwnHFDMzKwnHFDMzKwnHFDMzKwnajXKS9J+YGaWbEtJJk/WTV2vG3ztdbz2ul43zO/axyJi1pnhtQoo3ZA03c3wuKqp63WDr72O117X64b+XrubvMzMrCccUMzMrCccUE62Ne8C5KSu1w2+9jqq63VDH6/dfShmZtYTrqGYmVlP1DKgSFol6TlJuyTdmHH+VEl3p+cfk7Ri8KXsjy6u/ROSnpH0lKSHJM2yv255zHbtTfkulxSSKjEKqJvrlvTb6f/7Tkl/Mugy9ksXn/flkh6R9NfpZ/6SPMrZa5K2SXpZ0tNtzkvSH6f/Lk9Jek9PfnFE1OoBDAPfB34BWAR8Bzi/Jc8k8IX0+RXA3XmXe4DX/j5gcfr8ujpde5rvDOBbwKPAeN7lHtD/+Urgr4El6fHb8i73AK99K3Bd+vx8YE/e5e7Rtf8T4D3A023OXwJ8g2QH3PcCj/Xi99axhnIhsCsidkfEYeAuYHVLntXA9vT5vcBFkrK2Hi6bWa89Ih6JiMbG64+S7IhZBd38vwN8GvhPwN8OsnB91M11/0vg8xHxKkBEvDzgMvZLN9cewM+lz99Ch91fyyQivgUc7JBlNfCVSDxKsgvuWQv9vXUMKGcDLzYd703TMvNEstHXj0k28iq7bq692bUkdzFVMOu1S3o3cG5E/LdBFqzPuvk//yXglyT9b0mPSlo1sNL1VzfX/h+ANemGf18H/vVgipa7uf4t6MrI7FkqJ6um0TrUrZs8ZdT1dUlaA4wD/7SvJRqcjtcuaQi4BbhmUAUakG7+z0dImr1+k6RG+j8lvSsiftTnsvVbN9d+JXB7RPyRpF8H7kiv/Vj/i5ervvyNq2MNZS9wbtPxOZxczf27PJJGSKrCnaqPZdHNtSPpYmADcGlEvD6gsvXbbNd+BvAu4JuS9pC0K++oQMd8t5/3+yLijYj4AfAcSYApu26u/VrgHoCI+DbwJpK1rqquq78Fc1XHgPI4sFLSeZIWkXS672jJswNYmz6/HHg40p6skpv12tNmn1tJgklV2tJhlmuPiB9HxNKIWBERK0j6jy6NiLLvGd3N5/2/kgzGQNJSkiaw3QMtZX90c+0vABcBSHonSUDZP9BS5mMHcHU62uu9wI8j4qWFvmntmrwi4oik64H7SUaBbIuInZJuBqYjYgfwZZKq7y6SmskV+ZW4d7q89j8ETge+lo5DeCEiLs2t0D3S5bVXTpfXfT/wQUnPAEeBT0bEgfxK3RtdXvu/Bb4o6fdImnyuqcLNo6SvkjRhLk37h24CTgGIiC+Q9BddAuwCDgHrevJ7K/BvZ2ZmBVDHJi8zM+sDBxQzM+sJBxQzM+sJBxQzM+sJBxQzM+sJBxRrEzY2AAAAy0lEQVQzM+sJBxQzM+sJBxSzHEn61XQ/ijdJOi3dj+RdeZfLbD48sdEsZ5L+gGTJjzcDeyPiP+ZcJLN5cUAxy1m6ztTjJHuw/OOIOJpzkczmxU1eZvk7k2T9tDNIaipmpeQailnOJO0g2U3wPOCsiLg+5yKZzUvtVhs2KxJJVwNHIuJPJA0DfyXp/RHxcN5lM5sr11DMzKwn3IdiZmY94YBiZmY94YBiZmY94YBiZmY94YBiZmY94YBiZmY94YBiZmY94YBiZmY98f8BHCZ6zS4JSLEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x159112053c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.random.rand(100) # uniformly distributed x\n",
    "y = np.ones(x.shape)\n",
    "z = np.cos(2*x/0.25)\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(x,z,'ro')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "信号如图10.6所示。在图中，我们看到数据显示了一个很好的模式。模式的可见性来自数据确实有一些模式或抽样足够好的。正如我们在上一个例子中所看见的，原始变差图是非常杂乱的，所以我们只绘制平均变差图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'foo' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-4eac593f7550>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mfoo\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'foo' is not defined"
     ]
    }
   ],
   "source": [
    "foo = krige.OK(x,y,z)\n",
    "DE,GE = foo.variogram(var_type='averaged', min_lag=0.025)\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(DE,GE, 'r--o')\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.ylim((0,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "变差函数结果图如图10.7所示。变差图显示了清晰的趋势。变异性随着滞后距离的增加而增加。这意味着在更大的滞后距离上有更高的变异性。即该采样足够好捕捉数据的变异性。\n",
    "\n",
    "变差图在0.3的滞后距离附近稳定。这个距离大约是相关长度(在这个滞后距离滞后数据显示没有相互关系)。\n",
    "\n",
    "让我们再生成一个信号。图10.8显示了该信号。随着x的增加，变量(z)呈上升趋势，因此信号是非平稳的。这个变量变化的也很快，而且任何行为都不能从这个趋势中看出。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHh1JREFUeJzt3X2s5Fd93/H3Z9fe0Nvy5N0lRbbvvY7qSDGIBnxjQVUFUqAsW8mbSg6xtTY2MbmFyigqqIqtjQoyXZFQtaRpcciWODbsBYeQFLYEZDU8iCrg1Bc5OLaRYTH2emUaX68JInKIsfn2j99ce3Z2Hs7M/B7O7zeflzS68/C7M+fcmfv7zjnf86CIwMzMLMWOpgtgZmbt4aBhZmbJHDTMzCyZg4aZmSVz0DAzs2QOGmZmlsxBw8zMkjlomJlZMgcNMzNLdlbTBSjbnj17YnV1telimJm1yte+9rXHImLvpOM6FzRWV1fZ3NxsuhhmZq0i6aGU49w9ZWZmyRw0zMwsmYOGmZklc9AwM7NkDhpmZpas0aAh6WZJj0q6Z8JxPyfpaUmX1VU2M7NabGzA6irs2FH83NhoukRjNd3SuAXYN+4ASTuB3wJur6NAZma12diA9XV46CGIKH6ur2cdOBoNGhHxZeDxCYe9A/hj4NHqS2RmVqNDh+CJJ06/74knivsz1XRLYyxJ5wL/GvhQ02UxMyvdiRPT3Z+BrIMG8NvAr0fE0+MOkrQuaVPS5tbWVk1FMzOb0/LydPdnIPegsQbcJulB4DLgJkm/OHhQRByJiLWIWNu7d+LSKWZmeTh8GJaWTr9vaam4P1NZrz0VERdsX5d0C/CZiPhUcyUyMyvRwYPFz0OHii6p5eUiYGzfn6FGg4akjwOvAfZIOgm8GzgbICKcxzCz7jt4MOsgMajRoBERV0xx7DUVFsXMzBLkntMwM7OMOGiYmVkyBw0zM0vmoGFmZskcNMzMLJmDhpmZJXPQMDOzZA4aZmaWzEHDzMySOWiYmVkyBw0zM0vmoGFmZskcNMzMqrKxAaursGNH8TPjvb9TZb2fhplZa21swPr6s3uAP/RQcRtatRT6ILc0zMyqcOjQswFj2xNPFPe3WKNBQ9LNkh6VdM+Ixw9Kurt3+Yqkf1p3Gc3MZnLixHT3t0TTLY1bgH1jHv8O8OqIeBnwXuBIHYUyM5vb8vJ097dEo0EjIr4MPD7m8a9ExPd6N+8AzqulYGZm8zp8GJaWTr9vaam4v8WabmlM41rgc00XwswsycGDcOQIrKyAVPw8cqTVSXBoyegpSb9AETT++YjH14F1gOWWN/3MrEMOHmx9kBiUfUtD0suADwMHIuLUsGMi4khErEXE2t69e+stoJnZAsk6aEhaBv4EuCoivtl0eczMFl2j3VOSPg68Btgj6STwbuBsgIj4EPAfgN3ATZIAnoqItWZKa2ZmjQaNiLhiwuNvBd5aU3HMzGyCrLunzMwsLw4aZmaWzEHDzMySOWiYmVkyBw0zM+jk3hdVaMWMcDOzSnV074squKVhZtbRvS+q4KBhZlb33hct7gpz0DAzq3Pvi+2usIcegohnu8JaEjgcNMzM6tz7ouVdYQ4aZmZ17n3R8m1gPXrKzAzq2/tiebnokhp2fwu4pWFm3dCW5HLLt4F10DCz9mtTcrnl28AqIpouQ6nW1tZic3Oz6WKYWZ1WV4d3+ayswIMP1l2aVpL0tZT9itzSMLP2a3lyeWLX2sYG7NlTtEyk4npDrahGg4akmyU9KumeEY9L0u9IOi7pbkmvqLuMZtYCdc6zKNukrrWNDXjLW+DUqWd/59Qp+JVfaSRwNN3SuAXYN+bxNwIX9i7rwO/WUCYza5u6k8tlJt0nzds4dAh+9KMzf+/JJxuZ29Fo0IiILwOPjznkAPCRKNwBvEDSi+spnZm1Rp3J5bKT7pO61sZ1sTXQ/dZ0S2OSc4GH+26f7N1nZna6gweLpPePf1z8rGo0Utkzuid1rY3rYmug+y33oKEh950x3EvSuqRNSZtbW1s1FMusJdoyd6FNyki6978vf/u3sGvX6Y/3d60dPgxnn33mc+za1cjcjtyDxkng/L7b5wGPDB4UEUciYi0i1vbu3Vtb4cyy1qa5C3WbJ5jOm3QffF9OnSp+7t49vGvt4EH4gz8oHt+2ezfcfHMzczsiotELsArcM+KxfwV8jqLF8Urg/056vosvvjjMLCJWViKK09Hpl5WVpktWrqNHizpJxc+jRycfv7R0+t9kaWny75X1+/O8L9PWdQrAZqScs1MOquoCfBz4LvAjilbFtcDbgLf1HhfwQeDbwF8Ba5Oe00HDrEcafnKSmi5ZeWY5gZcRTOc5ec/6vswbrCZIDRqeEW7WVYswS3qWOu7YUZxyB0lFEr1qs74vFb+fnhFutuhavjBeklmS0k1PBJz1fclk1ruDhllXtXxhvLG2E9mjekrGBYCmg+ms70vTwW5bSh9Wmy7OaZh13LC+/Wn7+StMKFfGOY1qOKdh1nGj+vah+NZ++HA3WlPDbGwUkwhPnChaGCXWNTWn4aBhZu3SdCK7o5wIN7P2SZl0l0vf/oJy0DCzPKTOYG86kb3gHDTMLA+pCwF2eVRYCzinYWZ5cK6iUc5pmFm7NJGr8CrAU3PQMLM8NLH73mAO5corz9x/24HlNGc1XQAzM+DZnERF8xDOMCyHAsVS5evrz95eX3/2uO3kfH95F4xzGma2mEblULatrBQ/u77oY49zGmZm/Qa7mc45Z/zxJ06kLRK4YN1XDhpmXbNgJ7Ekw/IXP/jB8G1Uty0vT07OL+DuiI0GDUn7JN0v6bik64c8vizpi5LuknS3pP1NlNOsNRbwJJZkWP7iySfhec87fRvVbdsJ+EnJ+dS5JR3SWNCQtJNiV743AhcBV0i6aOCw3wA+EREvBy4Hbqq3lGYts4AnsSSjupkefxweewyOHh0+WXDSRMJM9rioU5Ojpy4BjkfEAwCSbgMOAPf1HRPA83rXnw88UmsJzdpmAU9iSZaXhye0t7uZtgPEMOMem/S8HdRk99S5wMN9t0/27uv3HuBKSSeBzwLvqKdoZi3lxfyGq2oOyAKug9Vk0NCQ+wbHv10B3BIR5wH7gY9KOqPMktYlbUra3NraqqCoZi2xgCexJFWtV7WA62A1Nk9D0quA90TEG3q3bwCIiPf1HXMvsC8iHu7dfgB4ZUQ8Oup5PU/DFl6FG/VYd7VhnsadwIWSLpC0iyLRfWzgmBPAawEk/QzwHMBNCbNxDh4sJp79+MfFz0UPGIswBLnGOjYWNCLiKeA64HbgGxSjpO6VdKOkS3uHvQv4VUlfBz4OXBNdm8K+SBbhn9fyMs8Q5LI/r1V9/useZp2ykXibLhdffPG0+6lbHY4ejVhaiig+1sVlaam436wqKyunf+a2Lysr439v3Of16NHi96XiZ8pnuMrP/6x1HABsRsI51mtPWT1WVxdmDZ9O6EpeZNY9OkZ9Xnfvhr/7u9PnwiwtTU5+V/n5L2kfkjbkNGyReP5Ae5TZ3dF0l+SsQ5BHfS5PnZpt8mSVn/+ah1k7aFg9PH+gPcqaVZ7DkiazDkGe9nM56eRf5ee/7mHWKX1Ybbo4p5Ep5zTaQxreRy5N9zwl9bXPrcwcxO7d5edIyjBLHQeQmNNo/CRf9sVBI2MlfLCtBmWd7MsKPk0Z9nmd5+Sf+effQcMsB5mfKIYq61txLi2NsrXxPU2QGjSc0zCrSg59+rMoa2mMri5psuCTJx00zKrS1DLlZYxYKuPEuIDrMi0CBw2zqjQxzDi31s128PnoR4vbV13l1QBazkHDrCpNDDPOcROm3ALZNJqeZ5IhBw2zqjTRp5/jJMocA1mKNge7CjlomFWliT79HCdR5hjIUrQ12FXMQcOsSgcPFi2L5eXiJHno0HTfVKftHslxxNKogLVjR97f2tsa7CrmoGFWpXmX5p72d3McsTQskAE8/XTe3T3jgt0i5zhSJnO06eLJfXaapidizTPBrUuT444ejdi5s131GTbJcfDSoaVw8NLotvC2v6lPu4x1meZZtrqkJa+z0cb69C8Rv2NH0Toa1JHl/VuxNLqkfZLul3Rc0vUjjnmTpPsk3SvpY3WX0Vosh0TmPInpHJPa85hUn7qHt6a8Xv8kx1GBbdFyHCnNkSouwE7g28BPAbuArwMXDRxzIXAX8MLe7RdNel53T9kzclgwb94F7rq0MvCw+px99rMrxw6+X2XXtb+rcvfuiF27pnu9LnUXDkHuCxYCrwJu77t9A3DDwDHvB946zfM6aNgzcvknnyevMnii27273QvlTTpxV/VepeQnJr1e14L4gFKDBvB5YP/AfUdSfnfMc14GfLjv9lXAfx845lO9wPHnwB3AvknP66Bhz+jSP3mX6rJtVFCvolWY8lopr9f0wIoKpQaN1JzGBcCvS3p3330TEyYTaMh9g1mysyi6qF4DXAF8WNILzngiaV3SpqTNra2tOYtlnZHj8NNZ5ZCfKVtKLqCs/E1q3mF5eXyuY96FHDuwLElq0Pgb4LXAT0r6X5KeX8JrnwTO77t9HvDIkGM+HRE/iojvAPdTBJHTRMSRiFiLiLW9e/eWUDTrjByXsZ7lxNHFiWaTAkKZkxJTgs/SEuzfX93SIV1ZliSlOQLc1Xf9GuCvgJMpvzvmOc8CHqBoxWwnwl8ycMw+4Nbe9T3Aw8Ducc/r7ilr3LgujFm7mXLJz5Rp2N9iOxledtfPuCR8//tU5d858/eQknMa/2bg9sXAzSm/O+F59wPfpBhFdah3343Apb3rAv4LcF8vUF0+6TkdNKxRk4LCrCeOLuY0IurNEaS8VpUj7nIYzTdGatDw5D6zMq2uFt0Og7YngM0zwa1/otnyctF1k0N3W5dMev9yfe4StGJyn1nnTMo9zDNhL8f8TNdUueBjjotJzsBBw6xMk4JC208cHRj984xhdUkZcTfr36Aro/lS+rDadHFOwxqVknto61j/LuVVZq1Ll/4GA8h9RnhVFwcNa1xbg0K/YXXIfPTPVGatS5f+BgNSg4a7pyxPbe4GmSb3kGM9R80nGJbEhXbOFZl13ksX58tMyUHD8tOVSVCT5FrPUbPPd+4cfnxKEj+34DjrgISurTw8i5TmSJsu7p7qgA53AZwmtZ51d3eNmk+w3X/fhTyAcxpnwDkNa63MJ0GVJqWedZ2k+gPTuB32ZglguX4JmDUYdyFnNYSDhrVXrieZsqXUs46/RdXbmtb5JaCjJ/Q6pAYN5zQsP22eyzBN331KPctMvI4q27AcBhQ5jDLmE9SVB8g1R9Q1KZGlTRe3NDqijd8YZ+lKmlTPsloa48pWdUugri62RWmhVgSvPWVWsyrWFtr+9tzfElhamv6b/7iyQfVrItWxbtY863qZ154yq10VY/jLWnpiXNnq6A6sY90sD4ethYOGWVmqOmmVccIdV7aurInU5lxYizhomJUl55PWpLJ1YQXdrgS/zDUaNCTtk3S/pOOSrh9z3GWSQtK8+5KbVSfnk1bOZStTF4Jf5hpLhEvaSbFr3+sp9gK/E7giIu4bOO65wJ9SbAl7XUSMzXI7EW6l8aZHtkDakAi/BDgeEQ9ExJPAbcCBIce9F3g/8MM6C9cKda/n0/96e/YUl1zWEiqbx/ybDdVk0DgXeLjv9snefc+Q9HLg/Ij4TJ0Fa4W6T2qDr3fqVHEp87VzWtRu1KJ9hw41Ux6zTDQZNDTkvmf6yiTtAD4AvGviE0nrkjYlbW5tbZVYxBKVfUKs+6T2a782fNZwWa+d2zf7KobP5hQUp9Xmslu5UmYAVnEBXgXc3nf7BuCGvtvPBx4DHuxdfgg8AqyNe94sZ4RXMSO27vV8xq1LVMZr5zabt+zy5Lw66qRZ6TmX3UpD7gsWAmcBDwAXUCS5vw68ZMzxX5oUMCLXoFHFCbHOk+yo1yrztXNb2bbsE2VuQXFbSj1zLbuVKjVoNNY9FRFPAdcBtwPfAD4REfdKulHSpU2VqxJVdHXUOSdg1I5tZb52brN5yx6imuuObyndnLmW3ZqRElnadGl1S2PaRfrqWNRv3IJ2/Zd5X7vrXSC5fltPaeHlWnYrFbl3T1V1yTJopJwQcz1ppnRNlXXyqCoI5rBibtve3/73NNeyW6kcNHJT1xLYZZvUysj95JHTCW/UZ6DJoJb698kh8FqlHDTaJrdE8LZxLY02nDxyC8aDJ9+3v735oOaAYOGg0T7Tntzq+kfP6Zv6LHIKxsP+lqPK13QL0xZOatDwKre5mGY0VJ0T4apa6K6uyWJ1bjU6qT7DRipFDH8+j0yyXKVEljZdWtvSiEhvPeTW5TKtOlsvdbxW6mukjEJr23tpnYG3e+2wtm9rWcW2qONUvVptan1GHSed/n7Osp2r2ZzasMqtzSq3iXCDJnXV1D1ZrOo9FlLrM6oL8m1v6/4+F9YZDhptlPMOcSn5ltyD3rRS6zMqP3TTTd44yFrDQaONct6FLWVZipyD3iymqY93lrOWc07DypWab+narnhdq48tnNSchoOGlavuJLeZlcKJcGtG17qezOw0Dho56NKuaDnnW3LVpfffOs/dU03bHm3Unzz2OP3F4fffMuGcRls4B7DY/P5bJlqR05C0T9L9ko5Lun7I4++UdJ+kuyV9XtJKE+WslHdFW2x+/61lGgsaknYCHwTeCFwEXCHpooHD7qLYF/xlwCeB99dbyhp0baJbLtqSJ/D7by3TZEvjEuB4RDwQEU8CtwEH+g+IiC9GxHZn7x3AeTWXsXoebVS+OlcBnpfff2uZJoPGucDDfbdP9u4b5Vrgc8MekLQuaVPS5tbWVolFrIFHG5UvZVZ6FWZp3fj9t5ZpLBEu6ZeAN0TEW3u3rwIuiYh3DDn2SuA64NUR8ffjnrd1ifBZeQbyaE2sAuxRUNZybUiEnwTO77t9HvDI4EGSXgccAi6dFDAWRpu6X5rQRJ6gqdaNWc2aDBp3AhdKukDSLuBy4Fj/AZJeDvweRcB4tNLStCVxCj5BTdJEnqBto6Da9Hm3vKTs1FTVBdgPfBP4NnCod9+NFEEC4M+Avwb+snc5Nuk5Z9q5r237YOe073Wu6tpDfVubdlNs2+fdaoF37ptC2yZYta28i6BNOQ1/fmyINuQ08tG2rgUP08xPm0ZBte3zbllx0ID2TbBq0wlqkbRlg6W2fd4tKw4a0M5v7m05QVl+2vh5t2w4aIC/udti8efd5uBEuJmZORFuM/L4fTMb46ymC2AZGRw2uj3THNx1YWaAWxrWzzPNzWwCBw17lsfvm9kEDhqLaljuwuP3zWwCB41FNGqV3P37PX7fzMZy0FhEo3IXn/2sx++b2VgOGm1R5lDYcbmLsmaae+iuWSc5aDRh2hNq2ZsunXPO8PvLyl14kyizznLQqNssJ9Qyh8JubMAPfnDm/WefPV/uoj8QXn21h+6adVSjy4hI2gf8V2An8OGI+M2Bx38C+AhwMXAK+OWIeHDcc2a/jMgsexmUuef1qNffvRsee2y659o2bC+JYarco9vM5pL9MiKSdgIfBN4IXARcIemigcOuBb4XEf8E+ADwW/WWsgKzzIWYNBR2mu6uUa/z+OOjf2eSYS2hYTx016z1muyeugQ4HhEPRMSTwG3AgYFjDgC39q5/EnitJNVYxvLNMhdi3FLW03Z3VTEXI2Xyn4fumnVCk0HjXODhvtsne/cNPSYingK+D+yupXSzSPnGP8teBuOWsp4231HFXgqjAs7OnR66a9Y1KRuJV3EBfokij7F9+yrgvw0ccy9wXt/tbwO7hzzXOrAJbC4vL5exx/r0jh6NWFqKKL7vF5elpeL+YceurERIxc9hx6SSTn/N7Ys0vqxlvf7286XW3cyyBGxGwrm7sUS4pFcB74mIN/Ru3wAQEe/rO+b23jFflXQW8P+AvTGm0I0lwmdJcLf5dQdtbBStmxMnipbH4cNuWZi1SPaJcOBO4EJJF0jaBVwOHBs45hhwde/6ZcAXxgWMRjW12F8uW3d6+1mzhdBY0IgiR3EdcDvwDeATEXGvpBslXdo77PeB3ZKOA+8Erm+mtAmaWuzPW3eaWY283WtZhs1VWFryCdzMWqEN3VPdsujf+L3WlNlC8HavZTp4cHGCRD9vE2u2MNzSsPl5m1izheGgYfPzNrFmC8NBw+bnbWLNFoaDhs0vl7kiZlY5Bw2b36KPHDNbIB49ZeVY1JFjZgvGLY1ced6DmWXILY0ced6DmWXKLY0ced6DmWXKQSNHnvdgZply0MiR5z2YWaYcNHLkeQ9mlikHjRyNmvcAHlFlZo1y0MjV4E54UIygeuihYhfu7RFVDhxmVqNGgoakcyT9b0nf6v184ZBjflbSVyXdK+luSb/cRFmz4RFVZpaBploa1wOfj4gLgc8zfBvXJ4A3R8RLgH3Ab0t6QY1lzItHVJlZBpoKGgeAW3vXbwV+cfCAiPhmRHyrd/0R4FFgb2Ulyn0GtkdUmVkGmgoaPxkR3wXo/XzRuIMlXQLsAr494vF1SZuSNre2tqYvzfYM7JzzBR5RZWYZUERU88TSnwH/eMhDh4BbI+IFfcd+LyLOyGv0Hnsx8CXg6oi4Y9Lrrq2txebm5nSFXV0tAsWglZVnk9A52NgochgnThQtjMOHvayImZVC0tciYm3ScZWtPRURrxv1mKS/lvTiiPhuLyg8OuK45wF/CvxGSsCYWVvyBV5J1swa1lT31DHg6t71q4FPDx4gaRfwP4GPRMQfVVoa5wvMzJI0FTR+E3i9pG8Br+/dRtKapA/3jnkT8PPANZL+snf52UpK43yBmVmSynIaTZkppwHOF5jZQms8p9E6zheYmU3kZUTMzCyZg4aZmSVz0DAzs2QOGmZmlsxBw8zMknVuyK2kLWDImiBJ9gCPlVicNnCdF8ci1nsR6wyz1XslIiYuCtu5oDEPSZsp45S7xHVeHItY70WsM1Rbb3dPmZlZMgcNMzNL5qBxuiNNF6ABrvPiWMR6L2KdocJ6O6dhZmbJ3NIwM7NkCxk0JO2TdL+k45KuH/L4T0j6w97jfyFptf5Sliuhzu+UdJ+kuyV9XtJKE+Us06Q69x13maSQ1IlRNin1lvSm3vt9r6SP1V3GsiV8vpclfVHSXb3P+P4mylkmSTdLelTSPSMel6Tf6f1N7pb0ilJeOCIW6gLspNhr/Kco9h3/OnDRwDH/FvhQ7/rlwB82Xe4a6vwLwFLv+tsXoc69454LfBm4A1hrutw1vdcXAncBL+zdflHT5a6hzkeAt/euXwQ82HS5S6j3zwOvAO4Z8fh+4HOAgFcCf1HG6y5iS+MS4HhEPBARTwK3AQcGjjkA3Nq7/kngtZJUYxnLNrHOEfHFiHiid/MO4Lyay1i2lPcZ4L3A+4Ef1lm4CqXU+1eBD0bE9wAiYuh2yy2SUucAnte7/nzgkRrLV4mI+DLw+JhDDlDsfBpRbJf9gt722nNZxKBxLvBw3+2TvfuGHhMRTwHfB3bXUrpqpNS537UU31DabGKdJb0cOD8iPlNnwSqW8l7/NPDTkv5c0h2S9tVWumqk1Pk9wJWSTgKfBd5RT9EaNe3/fZJF3IRpWIthcAhZyjFtklwfSVcCa8CrKy1R9cbWWdIO4APANXUVqCYp7/VZFF1Ur6FoUf4fSS+NiL+puGxVSanzFcAtEfGfJb0K+Givzj+uvniNqeQ8togtjZPA+X23z+PMpuozx0g6i6I5O64ZmLuUOiPpdcAh4NKI+PuaylaVSXV+LvBS4EuSHqTo8z3WgWR46uf70xHxo4j4DnA/RRBpq5Q6Xwt8AiAivgo8h2J9pi5L+r+f1iIGjTuBCyVdIGkXRaL72MAxx4Cre9cvA74QvcxSS02sc6+r5vcoAkbb+7hhQp0j4vsRsSciViNilSKPc2lEzLDBfFZSPt+fohj4gKQ9FN1VD9RaynKl1PkE8FoAST9DETS2ai1l/Y4Bb+6Nonol8P2I+O68T7pw3VMR8ZSk64DbKUZd3BwR90q6EdiMiGPA71M0X49TtDAub67E80us838C/hHwR72c/4mIuLSxQs8psc6dk1jv24F/Kek+4Gng30fEqeZKPZ/EOr8L+B+S/h1FF801Lf8iiKSPU3Qx7unlat4NnA0QER+iyN3sB44DTwBvKeV1W/53MzOzGi1i95SZmc3IQcPMzJI5aJiZWTIHDTMzS+agYWZmyRw0zMwsmYOGmZklc9Awq5ikn+vtZ/AcSf+wt4fFS5sul9ksPLnPrAaS/iPF0hX/ADgZEe9ruEhmM3HQMKtBb02kOyn27fhnEfF0w0Uym4m7p8zqcQ7F2l7PpWhxmLWSWxpmNZB0jGJHuQuAF0fEdQ0XyWwmC7fKrVndJL0ZeCoiPiZpJ/AVSf8iIr7QdNnMpuWWhpmZJXNOw8zMkjlomJlZMgcNMzNL5qBhZmbJHDTMzCyZg4aZmSVz0DAzs2QOGmZmluz/Azx7wIaq82UPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x159111c7898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.random.rand(100) # uniformly distributed x\n",
    "y = np.ones(x.shape)\n",
    "z = x + 0.2*np.random.randn(len(x))\n",
    "plt.clf()\n",
    "\n",
    "plt.plot(x,z, 'ro')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图10.8:z随着x的变化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foo = krige.OK(x,y,z)\n",
    "DE, GE = foo.variogram(var_type='averaged', min_lag=0.05)\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(DE,GE, 'r--o')\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.ylim(ymin=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "让我们现在估计它的变差函数。变差函数如图10.9所示。该变差函数有以下两个明显的特征:\n",
    "\n",
    "- 变差函数在0滞后距离时并不接近于零。\n",
    "- 变差函数随着滞后距离恒定增长。\n",
    "- 变差函数并不稳定。\n",
    "\n",
    "以上三个要素有如下揭示:\n",
    "\n",
    "- 数据的采样比数据的变异性大。\n",
    "- 数据的范围不足以捕捉数据的变异性。\n",
    "- 数据是非固定的。\n",
    "\n",
    "我们移除z中的趋势，然后绘制变差函数图。去趋势后的数据如图10.10所示，并且变差函数如图10.11所示。现在的数据没有趋势，因为趋势被移除了。而且该数据看上去只是没有趋势的随机数。变差函数结果图也表明，数据在所有的滞后距离处具有相同的变异性，这意味着它是随机行为。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAGnZJREFUeJzt3X+sZGV9x/HPdxdXvRV/sLtawo97McGmSBp1b420abVVK/LH0j9Mi7koENNtltIm1ZhCNqkGu8FojE2jCd5aUmSvFTBp3RoMVZTatEK5BKWi2bJSdtlKy4JKbDbKr2//OHNhmDs/npk55/lxzvuVTObO3HPnPmfmzPme5/l+z3PM3QUAQIgtqRsAACgHQQMAEIygAQAIRtAAAAQjaAAAghE0AADBCBoAgGAEDQBAMIIGACDYSakbULcdO3b40tJS6mYAQFHuvvvuR91956TlWhc0lpaWtL6+nroZAFAUMzsSshzDUwCAYAQNAEAwggYAIBhBAwAQjKABAAhG0ACAlNbWpKUlacuW6n5tLXWLxmpdyS0AFGNtTdqzRzpxonp85Ej1WJJWVtK1awx6GoinsCMqoHH79j0XMDacOFE9nyl6GoijwCMqoHFHj073fAboaSCOAo+ogMadeeZ0z2eAoIE4CjyiAhq3f7+0sPD85xYWquczRdBAHAUeUQGNW1mRVlelxUXJrLpfXc16yJaggTgKPKIColhZkR58UHrmmeo+44AhETQQS4FHVAA2o3oK8aysECSAwtHTAAAEI2gAAIIRNAAAwZIGDTM738wOmdlhM7tyzHLvMjM3s+WY7QMAPF+yoGFmWyV9WtI7JZ0j6d1mds6Q5U6W9CeS7ozbQgDAoJQ9jTdKOuzuD7j7E5K+IOnCIct9RNLHJP0sZuMAAJulDBqnSXqo7/Gx3nPPMrPXSzrD3b8cs2EAgOFSBg0b8pw/+0uzLZI+KekDE1/IbI+ZrZvZ+vHjx2tsIgCgX8qgcUzSGX2PT5f0w77HJ0s6V9LtZvagpDdJOjgsGe7uq+6+7O7LO3fubLDJANBtKYPGXZLONrOzzGybpIskHdz4pbs/7u473H3J3Zck3SFpt7uvp2kuACBZ0HD3pyRdIelWSd+XdJO732dmV5vZ7lTtAlqFqyWiZknnnnL3WyTdMvDcn49Y9i0x2gS0BldLRAM4IxxoK66WiAYQNIC24mqJaABBA2grrpaIBhA0gLbiaoloAEEDaCuulogGcOU+oM24WiJqRk8DABCMoAEACEbQAAAEI2gAQFNaOI0LiXAAaEJLp3GhpwEATWjpNC4EDQBoQkuncSFoAEATWjqNC0EDAJrQ0mlcCBqxtLCKAsAYLZ3GheqpGFpaRQFgghZO40JPI4aWVlEA6B6CRgwtraIA0D0EjRhaWkUBoHsIGjG0tIoCQPcQNGJoaRUFgO4haMSysiI9+KD0zDPVPQEDyAtl8UEouQUAyuKD0dMA0A7z9BQoiw9GTwNA+ebtKVAWH4yeRskYgwUq8/YUKIsPRtAo1caR1ZEjkvtzR1YEDnTRvD0FyuKDETRKxRhsGHpj3TBvT4Gy+GAEjZxMs4NjDHayrvbGuhgo6+gpUBYfxt1bddu1a5cX6cAB94UF92r3Vt0WFqrnh1lcfP6yG7fFxZitzlsX36Npt6M2OXCg+mzNqvsurHONJK17wD7WqmXbY3l52dfX11M3Y3pLS9WR8KDFxeqoZ9BgtYhUHVnRpX7Oli3VbnOQWXU02UbTbkdAj5nd7e7Lk5ZjeCoX0w43MQY7WRcrYhi2RMMIGrmYZQfHGOx4XayI6WKgRFRJg4aZnW9mh8zssJldOeT37zez75nZvWZ2m5ktpmhnFF3cwTWti70xtiM0LFnQMLOtkj4t6Z2SzpH0bjM7Z2CxeyQtu/uvSPqipI/FbWVEXdzBxZC6Nxa7kontCA1Llgg3s/Mkfdjd39F7fJUkufs1I5Z/vaRPufuvj3vdYhPhaB+KFVCQEhLhp0l6qO/xsd5zo7xP0lcabRFQJ07ARAulnLDQhjw3tNtjZhdLWpb05hG/3yNpjySdScIPuaCSCS2UsqdxTNIZfY9Pl/TDwYXM7G2S9kna7e4/H/ZC7r7q7svuvrxz585GGgtMjUomtFDKoHGXpLPN7Cwz2ybpIkkH+xfo5TE+oypgPJKgjcDsqGRCCyULGu7+lKQrJN0q6fuSbnL3+8zsajPb3Vvs45JeIulmM/u2mR0c8XJomzbMn0QlUxxt2FYKwjQiyA9VRwjFtlKbEqqnMI82H11RdYRQbCvRcbnXEs17acvcUXWEUGwr0dHTKFGqo6tYvRuqjhCKbSU6gsaGkoZ7UhxdxbygUduqjkratkqx8Z4eOVIVGfQrYVspeZsIuehGSbeZLsJU2oVrUlxcaNT/3Pi/db9XbbmgTmnbVmohn/uw99SsuW2xbpluEwq8CFPynXzdt5mCRmlXeEux0W18KUfdMtjos1TatpVS6HZdx3ua8qAk020iNGhQciuVeYW3tbUqh3H0aDV+u39/s0nwHTukxx4bvwxXh9usxG0rldCrDs77nqYu0810m6DkdholJtNiTvm9tib99KeTl6NiZbMSt61UQnN1876nqct0C98mCBpS+xKvddu3T3riicnLFbLRR8W2FS50Zzrve5q6TLf0bSJkDKuk20w5Dfe0Y5y5J30n5TPIaYyX++ebi2lydfO8pznkRDLcJkQivBCZVlI8z/bt4wPG9u15tRflirEznfc7V8J3dgahQYNEeGqhyb9U1takyy6Tnnxy9DK5tBUINU8hSe7f2RmFJsIJGqllWknxrFFfkH65tBWIIffv7IyonipF7pUUIcnBXNoKxJD7d7ZhBI3Ucq+kmPRFyKmtQAy5f2cbRtBILfcL9Qz7gmzM9ZNbW9toljmK6pjXqOS5kZqW+3e2aSHZ8pJuxVVPlSDD8sBOmKVKp47KnpZWB2E8UT0FFG6WKp06KntaWh2E8UiEA6nUNbQzy5nLdZztnPqMaWSNoAHUqc7rjsxSpVNHZU/Hq4OSKCiHRNBAPgr64oxU52R4s1Tp1FHZ0/HqoOhiXuCsDiGJj5JuJMIL1Zbk66h5usxme71ZihDqKFyg+CGeTK6vIRLhKEpbkq9tWQ/Ek8kZ5iTCUZa2JF9DhnbaMAyH+hSWQyJoIA+FfXFGmnTiV2nj12heYTkkggbyUNgXZ6xxV1VMfdU45KG/t7lvn3TJJePPMM+pdxqS+CjpRiK8YF1IvtadKEd5pi36GLa8mfvevbU2S4GJcHoayEfM656nMs8wXE5Hm5jdtL3NYcu7S9dem2QbIGgAMc06DEcupD2mLfoY9bx7kmFNgkapOOos06wzpJaaC2E73Wza3ua4XmiK6sKQMaySbp3IabTlRDiEKzEXUsJ2miKPNktOY9TnX+MJgArMaSTfydd960TQyOQMUkRU4meee5tTBrVpg9XevZsDR81tDQ0anBFeokzOIEVEGzmN/iGqhYW8L/6T+3Za2tn7a2vVcOTRo9WQ1f79tX72nBHeZm05EQ7hSrxaXG7b6WB+ZVjAkKrnc8y9ZFJdmDRomNn5ZnbIzA6b2ZVDfv9CM7ux9/s7zWwpfiszlPuJcCQ/m5HJTiNYTtvpsOqzjcsWD0Nl2mghY1hN3CRtlfQDSa+WtE3SdySdM7DM5ZKu7f18kaQbJ71uJ3Ia7vmeCBc6Tpxr+1GvXD7nUfmVUQnmnHIvkajOnIaZ3SbpE+5+S99zq+6+Z9ZgZWbnSfqwu7+j9/iqXhC7pm+ZW3vLfMvMTpL0P5J2+phGdyKnkbOQceISx+dRtlH5lXFyyb1EUndO4yxJf2ZmH+p7buKLT3CapIf6Hh/rPTd0GXd/StLjkrYPvpCZ7TGzdTNbP378+JzNwlxCTlAq9ZwDlGtUHmVxsbpN8zcdFxo0fiLprZJeZWb/aGYvq+F/DxtQHDwUCFlG7r7q7svuvrxz584amoaZhSQ/2zINOsoxLr8yKfeSU44uh7aEjGFJuqfv50sl/YekYyF/O+Y1z5N0a9/jqyRdNbDMrZLO6/18kqRHpWpIbdStMzmNXIXkNHKv348plzH/Lhj3Xo/6XU4nKDbcFtV5cp+kPxx4vEvSdSF/O+Y1T5L0gKqhr41E+GsHlvkjPT8RftOk1yVoZGDSjjCnL2LdpgkCbX4f2iKnA5yG21Jr0GjqJukCSf+pqopqX++5qyXt7v38Ikk3Szos6d8lvXrSaxI0CtHGI+xpg0BOOyQMl9P0LQ23JTRocEY4UJdpzzDO/Yxp5HXWeMNt4YxwYF7TJh2nTfDndsY0NsvpBMVM2kLQAIaZ5foV0waBTHYCUeRQ9TOLnKZvyaUtIWNYJd3IaaAWs+QbZklstzG3Myh1wr8L73ENVEIivIkbQQObzLLTmDXpyA5qs5QJ/9QBqyChQYNEONpt1ilLckqAli5lwp/PMRiJcECafcqSnPINpeYDNqRM+MeYfaD0z2dKBI0UOraRJTXrTiOXpOMsCfncpAzATQesNnw+0woZwyrpln1OgzHWuEo/ga709m9Iletp+vvWls/Hw3Ma9DTqFNKDYIbXuHIaZppFWyZ33LiA1A03VI/f8544veyme4xt+XymERJZSrol62mEHtHkNC1BV5Rc0dSiI9lW9rJb9PmInkZkoT0IzgKOr7TLpPYrvafUr4297DZ9PoEIGnUJ7aZ2cCPDHHJJyNehjUM5bfp8AnGeRl2mqQdfW6uOro4erXoY+/e3eiMDJHHOROY4TyO2aXoQswyXUKaL0tHLbgWCRl2a7KZ2sRYc7VPaUA4HakMxPFUCuvVAXLNOP1MwhqfapI0JRCBnbaz0qglBowSU6QJxcaA2EkGjBCQQgbg4UBuJoFGCmAlEkn/N4b0dL6f3hwO10UJOGy/plv2EhTlr4zQPueC9HW5jipeNqXRyen9GTT9T8rQ0Y4gr92FqLZpHJzulvbcxdpjDAmnu70+Lgz9BA+MN+/IzmWJzSnpvR+0Y9+6td4c5KpDm/P6UFvynEBo0OE+ji0bVoL/4xdJjj21envNB5lfSuTaj2rp1q/T005ufn3UdRl0Gto7XbkrKS9c2jPM0MNqoGnSJ5F9TckishiaaR5WVDgsY45afZFIlUo7bHlVVBI1WmLbqZNSX/Ec/Kmuah5KknkJj0lQ0/dvQlhG7ha1bhz8/6w5zWCA1q+5z3fZyCP6phYxhlXTrXE5jlsRci8dlMcK4zzwkId1ETsO9zEqkEtscQCTCO2KWANDiChCMMC4RP2ob2rq1M+WmUWT+3oUGDRLhpZs1Mcc1PbplXCL+6NHWJnezUcAEiCTCu2KWxBwBo345nc08zLixeJK7zWvTBIgh3ZGSbp0bnpp2qImhqfqV8p6OO2GvhPaXrIDzdEROo0OmGSvtShI85vjxNO9pruPaubarLQr43oUGDXIaXdPik5OeFXv8OPQ9LWBcGw0p4LPPOqdhZqeY2VfN7P7e/SuGLPM6M/uWmd1nZvea2e+naGvrdGH8Ovb4ceh72qZxbUwn9Xk6NUqVCL9S0m3ufrak23qPB52Q9F53f62k8yX9pZm9PGIb26kLJyfFvoBO6Hs66v8fOZJ3Er3L6ixwWFmppkR55pnqvsCAISlNTkPSIUmn9n4+VdKhgL/5jqSzJy3XyZzGtNo+fp1i/DjkPR3VrtymBEelYwUCyjmnYWY/cfeX9z3+sbtvGqLq+/0bJV0v6bXuPnbgnZwGsh0/HtYus+H5kNwm6uuikiaZrEHynIaZfc3MvjvkduGUr3OqpBskXTYqYJjZHjNbN7P148ePz9bg3OvsES7X8eNh7Rp10Ma1qNPjOuFDpeppHJL0Fnd/uBcUbnf3Xxqy3Esl3S7pGne/OeS1Z+pp5Hpkivbr2NFsUTr22STvaUxwUNIlvZ8vkfSlwQXMbJukv5f0udCAMTOqWpBKFwoTpDx68oNtuPzy8W3qymczrZDER903SdtVVU3d37s/pff8sqTP9n6+WNKTkr7dd3vdpNeeKRFewNmaaLG2FybkkFAOncl3sE1t/2z6KOdEeJNmGp7qWDcUiCqH79eoNgzq8Hc+9+GpvNANBZqTQ0I59H+VmuSOOPxH0JDyrbYB2iCHWQhC/9ek5XLIzQyadFXGmhE0NrTlbM1c5PjlQho59OSHtWHQpDZF3jkHi13IE5L4KOnGGeGJ9CcMt293f8EL0iY+kZccEsqDbdi7d7o25TpTbU2FPCIRXriSLpQ07DyXYTqcZEQL5DpDdE2FBiTCS5ZrN3iUYd3jYUpNMrYdQ4lhcsjNDBN5+I+gkaPSTjYMDQapv1zYrLQDlJRyyM0ME7mQh+GpHOXaDR4lpAaeaVnylMM5FCUpadh4SgxPlSzXbvAow47Atm2Ttm+nhDl3OZxDkZtxw3VUWRI0spRrN3iUYd3j666THn2001+uIpR2gNI0husmImjkqMSTDTkCK1PoAUpXkuWl5RMTIGjkip1wd6TcIYccoHTp6JvhuolIhAMplXAtly4ly7u0rgNIhAOxzdJjKGE4pEtH36XlExMgaAB1mHUIp4QdcluT5cOCfIn5xMgYngLqMOuwRgnDISUMoU2rjes0J4angJhm7TGUMBzSxqPvEoYFM0XQAOow6xBOKTvktlXzlTAsmCmCRkm6Uitfonl6DG3bIZegrXmaCAgapehSrXyJSukxoFLCsGCmCBqlyH0Mll4QPYYchG6HBPmZUT1VipxnvqUSBTlgO5xLaPUUQaMUOZdm5tw2dAfb4VwouW2bnMdgqURBDtgOoyBolCLnMVgqUZADtsMoCBolyTXRmnMvCN0Razucp+ijDQUj7t6q265duxwJHDjgvrjoblbdHziQukXooqa3wwMH3BcW3KuylOq2sBD2f+b52wgkrXvAPpZEOACEmifZnnminkQ4ANRtnmR7SxL1BA0ACDVPsr0liXqCBgCEmifZ3pKCEYIGAISap/Q957L5KZAIBwCQCAcA1C9J0DCzU8zsq2Z2f+/+FWOWfamZ/beZfSpmGwEAm6XqaVwp6TZ3P1vSbb3Ho3xE0j9HaRUAYKxUQeNCSdf3fr5e0u8OW8jMdkl6laR/itQuAMAYqYLGq9z9YUnq3b9ycAEz2yLpE5I+GLltAIARTmrqhc3sa5J+ccivQi81d7mkW9z9ITOb9L/2SNojSWcWdqIMAJSksaDh7m8b9Tsz+18zO9XdHzazUyU9MmSx8yT9hpldLuklkraZ2f+5+6b8h7uvSlqVqpLbetYAADCosaAxwUFJl0j6aO/+S4MLuPuzZ7yY2aWSlocFDABAPKlyGh+V9HYzu1/S23uPZWbLZvbZRG0CAEzAGeEAAM4IBwDUj6ABAAhG0AAABCNoAACCETQAAMEIGkAXra1JS0vSli3V/dpa6hahEKlO7gOQytqatGePdOJE9fjIkeqxVNxV5BAfPQ2ga/btey5gbDhxonoemICgAXTN0aPTPQ/0IWgAXTNqJmhmiEYAggbQNfv3SwsLz39uYaF6HpiAoAF0zcqKtLoqLS5KZtX96ipJcAShegroopUVggRmQk8DABCMoAEACEbQAAAEI2gAAIIRNAAAwVp3uVczOy7pyIx/vkPSozU2pwSsc3d0cb27uM7SbOu96O47Jy3UuqAxDzNbD7lGbpuwzt3RxfXu4jpLza43w1MAgGAEDQBAMILG862mbkACrHN3dHG9u7jOUoPrTU4DABCMngYAIFgng4aZnW9mh8zssJldOeT3LzSzG3u/v9PMluK3sl4B6/x+M/uemd1rZreZ2WKKdtZp0jr3LfcuM3Mza0WVTch6m9nv9T7v+8zs87HbWLeA7ftMM/uGmd3T28YvSNHOOpnZdWb2iJl9d8Tvzcz+qvee3Gtmb6jlH7t7p26Stkr6gaRXS9om6TuSzhlY5nJJ1/Z+vkjSjanbHWGdf0vSQu/nvV1Y595yJ0v6pqQ7JC2nbnekz/psSfdIekXv8StTtzvCOq9K2tv7+RxJD6Zudw3r/ZuS3iDpuyN+f4Gkr0gySW+SdGcd/7eLPY03Sjrs7g+4+xOSviDpwoFlLpR0fe/nL0p6q5lZxDbWbeI6u/s33H3jwtF3SDo9chvrFvI5S9JHJH1M0s9iNq5BIev9B5I+7e4/liR3fyRyG+sWss4u6aW9n18m6YcR29cId/+mpB+NWeRCSZ/zyh2SXm5mp877f7sYNE6T9FDf42O954Yu4+5PSXpc0vYorWtGyDr3e5+qI5SSTVxnM3u9pDPc/csxG9awkM/6NZJeY2b/amZ3mNn50VrXjJB1/rCki83smKRbJP1xnKYlNe33PkgXL8I0rMcwWEIWskxJgtfHzC6WtCzpzY22qHlj19nMtkj6pKRLYzUokpDP+iRVQ1RvUdWj/BczO9fdf9Jw25oSss7vlvS37v4JMztP0g29dX6m+eYl08h+rIs9jWOSzuh7fLo2d1WfXcbMTlLVnR3XDcxdyDrLzN4maZ+k3e7+80hta8qkdT5Z0rmSbjezB1WN+R5sQTI8dPv+krs/6e7/JemQqiBSqpB1fp+kmyTJ3b8l6UWq5mdqs6Dv/bS6GDTuknS2mZ1lZttUJboPDixzUNIlvZ/fJenr3sssFWriOveGaj6jKmCUPsYtTVhnd3/c3Xe4+5K7L6nK4+x29/U0za1NyPb9D6oKH2RmO1QNVz0QtZX1Clnno5LeKklm9suqgsbxqK2M76Ck9/aqqN4k6XF3f3jeF+3c8JS7P2VmV0i6VVXVxXXufp+ZXS1p3d0PSvobVd3Xw6p6GBela/H8Atf545JeIunmXs7/qLvvTtboOQWuc+sErvetkn7HzL4n6WlJH3T3x9K1ej6B6/wBSX9tZn+qaojm0sIPBGVmf6dqiHFHL1fzIUkvkCR3v1ZV7uYCSYclnZB0WS3/t/D3DQAQUReHpwAAMyJoAACCETQAAMEIGgCAYAQNAEAwggYAIBhBAwAQjKABNMzMfrV3PYMXmdkv9K5hcW7qdgGz4OQ+IAIz+wtVU1e8WNIxd78mcZOAmRA0gAh6cyLdpeq6Hb/m7k8nbhIwE4angDhOUTW318mqehxAkehpABGY2UFVV5Q7S9Kp7n5F4iYBM+ncLLdAbGb2XklPufvnzWyrpH8zs99296+nbhswLXoaAIBg5DQAAMEIGgCAYAQNAEAwggYAIBhBAwAQjKABAAhG0AAABCNoAACC/T/AYJsy5NZ7wwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x159111c7dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = z-x\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(x,z, 'ro')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foo = krige.OK(x,y,z)\n",
    "DE, GE = foo.variogram(var_type='averaged', min_lag=0.05)\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(DE,GE, 'r--o')\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.ylim(ymin=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，让我们对站点数据执行克里金插值，来绘制变量地图和非确定性估计。首先我们估计经验变差函数，然后通过命中和试验，将理论变差函数(在这个情况下为spherical)拟合。经验和理论拟合函数如图10.12所示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([6.8, 4.3, 5.9, 11.6, 5.5, 10.8, 8.6, 12.6, 14.7, 13.9,\n",
    "9.47, 14.3, 8.9, 11.9, 11.75])\n",
    "y = np.array([6.4, 5.0, 6.0, 4.9, 2.7, 8.2, 3.9, 6.7, 10.4, 10.9, 5.6,\n",
    "11.0, 7.3, 6.7, 10.8])\n",
    "z = np.array([10, 11, 12, 9, 12, 8, 10, 8, 6, 6,\n",
    "10, 6, 8, 8, 6])\n",
    "foo = krige.OK(x,y,z)\n",
    "DE, GE = foo.variogram(var_type='averaged', min_lag=0.5)\n",
    "\n",
    "model_par = {'nugget':0, 'range':20, 'sill':6}\n",
    "\n",
    "lags = np.linspace(0,6)\n",
    "G = foo.vario_model(lags, model_par, 'spherical')\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(DE,GE, 'rs', ms=10)\n",
    "plt.plot(lags, G, 'g', lw=3)\n",
    "plt.xlabel('lag distance')\n",
    "plt.ylabel('variogram')\n",
    "plt.ylim(ymin=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，首先我们生成所需的网格，在其上插入数据。这不只是网格格式，它可以在任何位置。但是为了制图的目的，需要网格格式的位置。`OK`类中的`krige`方法用于执行克里金插值。`krige`方法分别为Krigged值和方差生成两个属性`Zg`和`s2_k`。这两个属性都是一维数组，所以为了制图，我们需要将他们的形状转换为二维数组。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xg, Yg = np.meshgrid(np.linspace(4,16), np.linspace(2,12))\n",
    "foo.krige(Xg, Yg, model_par, 'spherical')\n",
    "krig_z = foo.Zg\n",
    "var_z = foo.s2_k\n",
    "krig_z.shape = 50,50\n",
    "var_z.shape = 50,50\n",
    "\n",
    "plt.clf()\n",
    "plt.pcolor(Xg, Yg, krig_z)\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.clf()\n",
    "plt.pcolor(Xg, Yg, var_z)\n",
    "plt.plot(x,y, 'ro', ms=12)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图10.13显示了Krigged值的地图。图10.14显示了原始站点位置的变化情况。从图中可以很明显看出，站点附近误差的变化很小，而且它随着我们远离站点而增加，当周围都没有站点时，误差便得非常高。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
